基于環(huán)形掃面測量的三維直流電阻率法任意各向異性模型響應(yīng)特征

殷長春，楊志龍，劉云鶴，張 博，齊彥福，曹曉月，邱長凱，蔡 晶

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

0 引言

直流電阻率法已廣泛應(yīng)用于地球淺部礦產(chǎn)資源勘查、水文地質(zhì)調(diào)查和環(huán)境工程領(lǐng)域[1]。它利用地表觀測的電勢分布，推測地下結(jié)構(gòu)的電阻率分布[2]。隨著野外采集儀器的快速發(fā)展以及數(shù)值算法和理論的不斷完善，高精度和高效地求解復(fù)雜地電模型，尤其是任意各向異性地下模型的三維直流電阻率響應(yīng)特征成為研究熱點。目前正反演算法采用的電各項同性介質(zhì)模型造成反演結(jié)果與地下實際地質(zhì)情況存在很大差異[3]，尤其是遇到具有層理面或裂縫的巖石等有方向依賴性的物質(zhì)結(jié)構(gòu)時，直流電阻率模型必須考慮地下電阻率各向異性[4-6]。

目前，人們已獲得均勻半空間傾斜各向異性介質(zhì)的電阻率響應(yīng)解析解[7]，而層狀各向異性介質(zhì)的電阻率響應(yīng)已由Li等[8]和Yin等[4,9]給出。對于復(fù)雜的各向異性模型，目前的數(shù)值解法主要包括：積分方程法[10]、有限差分法[11]和有限元法[12]。其中：Li等[10]利用積分方程法計算了任意地下各向異性異常體的地面電勢響應(yīng)；王威等[12]利用有限元法在直流電阻率法正演模擬中計算各向同性介質(zhì)中存在各向異性異常體時的電阻率響應(yīng)。在有限元方法的應(yīng)用中，非結(jié)構(gòu)網(wǎng)格以其可以模擬任意形狀異常體的突出優(yōu)勢備受重視，本文選擇了基于非結(jié)構(gòu)網(wǎng)格的自適應(yīng)有限元算法來模擬地下任意各向異性異常體的電阻率分布特征。為保證源附近及電阻率差異較大區(qū)域的模擬精度，我們在非結(jié)構(gòu)網(wǎng)格的基礎(chǔ)上采用了基于梯度恢復(fù)的局部網(wǎng)格自適應(yīng)加密技術(shù)，這不僅解決了全局網(wǎng)格自適應(yīng)造成的內(nèi)存較大、計算緩慢等問題，還可以根據(jù)經(jīng)驗對不同計算區(qū)域設(shè)定不同的單元誤差閾值，提高了網(wǎng)格自適應(yīng)加密的靈活性。

傳統(tǒng)的電法勘探中采用的測量方式是只觀測電場的一個同向分量，方位上的信息(垂直測線)往往被忽略。除此之外，視電阻率各向異性反?，F(xiàn)象的存在使得單一測線所得到的視電阻率不能準(zhǔn)確反映該方向的電阻率信息[4]。對于這種測量上的困擾，Yin等[13]在一維各向異性層狀介質(zhì)的視電阻率計算中采用了環(huán)形掃面測量方式，通過視電阻率極性曲線準(zhǔn)確地反映了地下各向異性介質(zhì)的方位信息。本文將該測量方法應(yīng)用于三維各向異性地下介質(zhì)的正演模擬中。利用二極裝置，以點源為圓心，以收發(fā)距為半徑，環(huán)形掃面測量多組共收發(fā)距各方位的電位信號，并且通過改變收發(fā)距以獲得不同深度的地電信息，在實際應(yīng)用中有很強(qiáng)的可操作性和靈活性。本文以淺部地表為研究目標(biāo)，首先計算任意各向異性半空間模型電阻率響應(yīng)，研究主軸電阻率旋轉(zhuǎn)過程中視電阻率極性曲線隨各向異性的分布特征；進(jìn)而，我們在任意各向異性半空間中加入了任意各向異性異常體，通過對幾種典型復(fù)雜各向異性模型的視電阻率分布特征進(jìn)行分析，總結(jié)了識別圍巖和異常體各向異性主軸電阻率和主軸方向的方法技術(shù)。

1 正演理論

1.1 電導(dǎo)率張量

對于任意各向異性地下介質(zhì)，電導(dǎo)率σ(σ=ρ-1，ρ為電阻率張量)是一個含有6個獨立分量的對稱正定張量，即

(1)

式中：x、y、z分別代表直角坐標(biāo)系的3個方向；σxy=σyx，σxz=σzx，σyz=σzy。根據(jù)Yin[6]，該電導(dǎo)率張量可由主電導(dǎo)率張量經(jīng)過3次歐拉旋轉(zhuǎn)得到，即

σ=RTσcR。

(2)

式中：σc為主電導(dǎo)率張量，

(3)

(4)

式中，α、β、γ分別是電導(dǎo)率主軸繞x、y、z軸的旋轉(zhuǎn)角[6]。

1.2 基于二次場的邊值問題

在三維半空間區(qū)域Ω中，地面Γ0上一點電源激發(fā)產(chǎn)生的電位u滿足Poisson方程[1]，地面邊界處的電位滿足Newman邊界條件，區(qū)域其他邊界Γ處的電位滿足第三類邊界條件，則二次場滿足的邊值條件表示如下：

·(σ·us)=-·((σ-σ0)·up)，r∈Ω。

(5)

(6)

σ·up·n-

(7)

式中：up和us分別表示背景場電位和由異常體激發(fā)的二次場電位；r表示空間任意一點指向發(fā)射源的距離向量；σ0表示圍巖電導(dǎo)率；n表示空間邊界的法向分量；B和BP是中間量，B=rT·σ-1·r，BP=rT·σ0-1·r。一次場的解析解為[8]

(8)

在Hilbert函數(shù)空間H1(Ω)內(nèi)取測試函數(shù)N，由Galerkin有限單元法可得

?Ω(·(σ·us)+

(9)

我們使用開源代碼TetGen[14]將區(qū)域剖分成Delaunay四面體網(wǎng)格，利用格林恒等式化簡式(9)后代入邊界條件，并對二次場在各單元中進(jìn)行插值[15]，得到線性方程組[16]

KUs=b。

(10)

式中：K是全局系數(shù)矩陣；b是右端源項；Us是待求的二次電位列向量。使用并行直接求解器MUMPS求解大型稀疏線性方程組可得到二次場電位，最后與背景場電位相加即可得到網(wǎng)格節(jié)點的總電位。在獲得測點的總電位后，利用文獻(xiàn)[17]給出的視電阻率定義即可計算相應(yīng)的視電阻率。

1.3 基于梯度恢復(fù)的誤差估計和局部網(wǎng)格自適應(yīng)加密策略

Zienkiewicz等[18]于1987年提出了基于梯度恢復(fù)的后驗誤差估計算法(又稱為Z-Z方法)。該方法實現(xiàn)起來簡單、高效，不依賴于特定問題，故算法的可移植性較高，與非結(jié)構(gòu)網(wǎng)格結(jié)合表現(xiàn)出很好的穩(wěn)健性(robustness)[19]。

在L2范數(shù)下，三維恒定電流有限元計算中的單元誤差可以表示為[17]

i=1, 2,,n。

(11)

式中：i表示單元編號；n表示目標(biāo)區(qū)域剖分的單元總數(shù)；Ωi表示單元子區(qū)域；uF為有限元解插值得到的單元電位梯度；du是利用超收斂路徑恢復(fù)技術(shù)[18]得到的與電位梯度精確值最為接近的單元梯度恢復(fù)值，即

du=Ga。

(12)

式中：

(13)

單元梯度恢復(fù)值在單元鄰域(element patch)內(nèi)求得，該鄰域由與單元Ωi相鄰的所有單元組成。通過求取式(13)中的a，可得到單元的恢復(fù)梯度值du，進(jìn)而可求得區(qū)域內(nèi)的單元誤差‖e‖i。

在自適應(yīng)網(wǎng)格加密策略中，我們首先設(shè)定每個目標(biāo)區(qū)域的最大誤差閾值:

‖e‖j；j=1,2，…,k。

(14)

式中，k表示目標(biāo)區(qū)域數(shù)量。進(jìn)而，可以求出每個單元一次迭代時的最大體積：

(15)

式中：Vi*表示單元原體積。式(15)表明自適應(yīng)網(wǎng)格的體積可通過單元誤差與最大誤差比值進(jìn)行調(diào)節(jié)。新的單元尺寸確定后，通過Delaunay網(wǎng)格剖分獲得新的迭代網(wǎng)格，然后重復(fù)上述步驟直到單元誤差小于對應(yīng)區(qū)域的單元誤差閾值為止。

由于全局自適應(yīng)加密方法計算的是整個模型內(nèi)的單元后驗誤差，在自適應(yīng)過程中會對整個模型區(qū)域進(jìn)行加密，而實際用到的數(shù)據(jù)只是整個區(qū)域中的部分區(qū)域，因此大量的網(wǎng)格數(shù)據(jù)不僅占用內(nèi)存，而且會極大增加計算時間。本文在此基礎(chǔ)上提出了局部自適應(yīng)加密，將模型區(qū)域分為計算區(qū)域和擴(kuò)邊區(qū)域。擴(kuò)邊區(qū)域體積較大，單元剖分對計算結(jié)果影響較小，所以對擴(kuò)邊區(qū)域利用單元誤差閾值直接代替單元誤差，這樣導(dǎo)致在網(wǎng)格的重新生成中該區(qū)域的單元基本保持原來的尺寸不變。對于計算區(qū)域，可以根據(jù)需要分為幾個不同的區(qū)域，根據(jù)經(jīng)驗對不同區(qū)域設(shè)定不同的單元誤差閾值，從而使得在網(wǎng)格重新生成的迭代過程中，不同區(qū)域的加密可不同，這樣不僅可以減少不必要的內(nèi)存浪費，還可有效地解決全局自適應(yīng)加密的盲目性。

2 精度驗證

本文算例運行平臺為Intel(R) Xeon(R) CPU E5-2667 v3 @ 3.20GHz 3.20 GHz。我們首先采用如圖1所示的任意各向異性的兩層模型來驗證本文算法的精度。為方便起見，采用電阻率進(jìn)行討論。第一層介質(zhì)的主軸電阻率及歐拉角為ρx1=400m，ρy1=100m，ρz1=100m，α1=0°，β1=0°，γ1=0°；第二層介質(zhì)的主軸電阻率及歐拉角為ρx2=40m，ρy2=10m，ρz2=10m，α2=0°，β2=0°，γ2=0°。假設(shè)點電源位于坐標(biāo)原點，采用二極裝置，通過與Wait[20]給出的半解析解對比，得到測量剖面的視電阻率和相對誤差曲線(圖2)。

圖1 兩層各向異性地電模型Fig.1 Illustration of the two-layer model

通過分析圖2可以得出如下結(jié)論：

1)有限元與半解析解的對比結(jié)果表明，本文算法具有較高的精度，能很好地滿足正演模擬需求。

2)收發(fā)距較小時，淺部信息得到較好的反映；隨著收發(fā)距增大，第二層電阻率得到越來越明顯的反映。

3)在收發(fā)距較小時，x、y方向的觀測結(jié)果發(fā)生分離，淺層的各向異性信息表現(xiàn)出來；而在收發(fā)距較大時，x、y方向的視電阻率也出現(xiàn)了分離，深部的各向異性信息表現(xiàn)出來。因此，通過改變收發(fā)距，層狀大地各向異性特征可以求解。

4)在收發(fā)距小時，x測線方向視電阻率趨近于沿層理的電阻率ρL1=100m，y方向趨近于沿層理電阻率ρL1和垂直層理電阻率ρT1的幾何平均值m；在收發(fā)距較大時，x測線方向視電阻率趨近于ρL2=10m，y方向視電阻率趨近于m。這實際上是典型的視電阻率各向異性反?，F(xiàn)象[10](anisotropy paradox)，即真實電阻率大的方向視電阻率小，而真實電阻率小的方向視電阻率大。

3 數(shù)值實驗

3.1 模型建立

下面我們重點研究任意各向異性圍巖和異常體視電阻率響應(yīng)特征。為此，設(shè)定如圖3所示的模型，建立z軸向下的直角坐標(biāo)系，坐標(biāo)原點位于地面中心，點電源位于坐標(biāo)原點處;異常體的長寬高均為100 m，埋深5 m，主軸電阻率為ρx1=10m，ρy1=40m，ρz1=10m；圍巖的主軸電阻率為ρx0=100m，ρy0=400m，ρz0=100m。在局部網(wǎng)格自適應(yīng)加密過程中，我們選出600 m600 m600 m的計算區(qū)域，設(shè)定單元相對誤差閾值為5%；

r. 收發(fā)距。a、b測線沿x軸；c、d測線沿y軸。圖2 本文算法結(jié)果與半解析解對比及相對誤差Fig.2 Comparisons of the result computed by this paper with semi-analytical solutions and the relative deviations

圖3 基于二極裝置的環(huán)形掃面測量與正演模型Fig.3 Concentric circular scanning measurement based on bipole configuration and 3D model

由于計算二次場，故對源的鄰域和異常體所在區(qū)域設(shè)定單元誤差為1%，計算區(qū)以外的擴(kuò)邊區(qū)域(范圍為2 000 m2 000 m2 000 m)不加密。由于計算區(qū)域相對于整個模型區(qū)域很小，為了便于查看網(wǎng)格加密情況，本文僅展示600 m600 m300 m的區(qū)域內(nèi)的y-z平面自適應(yīng)加密網(wǎng)格(圖4)。

對比圖4中a、b、c可知，在局部自適應(yīng)網(wǎng)格加密過程中，計算區(qū)域內(nèi)部網(wǎng)格不斷加密(而擴(kuò)邊區(qū)內(nèi)部的網(wǎng)格基本不變)，并且在第三次迭代后，計算區(qū)域的單元相對誤差和異常體內(nèi)部的單元相對誤差已經(jīng)達(dá)到要求。另外，觀察圖4c可以發(fā)現(xiàn)，由于設(shè)定的閾值各不相同，發(fā)射源附近區(qū)域和異常體內(nèi)部加密程度要大于計算區(qū)內(nèi)其他區(qū)域；異常體與圍巖的交界處由于電阻率差異較大導(dǎo)致單元相對誤差較大，網(wǎng)格加密程度較高。這是因為本文計算的是二次場，而二次場是由異常體與圍巖交界處電阻率的突變所激發(fā)，因此異常體的表面相當(dāng)于二次場的激發(fā)源，所以異常體表面單元相對誤差比較大，對應(yīng)的自適應(yīng)網(wǎng)格加密比較明顯。

a.初始網(wǎng)格及單元相對誤差；b.第一次自適應(yīng)迭代加密后的網(wǎng)格剖分和單元相對誤差；c.第三次自適應(yīng)迭代加密后的網(wǎng)格剖分和單元相對誤差，圖中白色虛線方框內(nèi)為異常體。圖4 三維模型局部自適應(yīng)加密網(wǎng)格Fig.4 Local adaptive mesh refinement for 3-D model

3.2 各向異性半空間模擬

我們首先研究各向異性半空間的視電阻率響應(yīng)特征。為此，我們在正演模擬中假設(shè)異常體和圍巖具有相同的各向異性，此時半空間的主軸電阻率為ρx/ρy/ρz=100/400/100m。圖5為半空間主軸電導(dǎo)率繞z軸和x軸旋轉(zhuǎn)時的視電阻率極性圖。其中，由原點到各測點的距離表示視電阻率大小，而對應(yīng)的方向表示實際測線方向。為簡化起見，針對不同各向異性的網(wǎng)格圖沒有展示。

由圖5a我們發(fā)現(xiàn)，視電阻率極性圖呈現(xiàn)橢圓狀，在主軸電阻率小的方向被拉伸(對應(yīng)橢圓長軸)，視電阻率大小為x與y方向的實際電阻率的幾何平均值；在電阻率大的方向被壓縮(對應(yīng)橢圓短軸)，視電阻率為x方向的實際電阻率。這正好與圖2中的視電阻率各向異性反?，F(xiàn)象相一致。利用這種視電阻率反?，F(xiàn)象可以有效地識別目標(biāo)體的各向異性特征。由圖5b可以發(fā)現(xiàn)，當(dāng)主軸電阻率繞x軸旋轉(zhuǎn)時，視電阻率極性圖也發(fā)生了改變：當(dāng)旋轉(zhuǎn)角度為0°時，地層層理直立，主軸電阻率為ρx/ρy/ρz=ρL/ρT/ρL=100/400/100m，x方向視電阻率等于m，y方向視電阻率等于ρx=100m；隨著層理的傾斜，橢圓型曲線的短軸方向逐漸被拉伸，當(dāng)旋轉(zhuǎn)至90°時，主軸電阻率為ρx/ρy/ρz=ρL/ρL/ρT=100/100/400m，地層變成水平各向同性，視電阻率曲線呈現(xiàn)出圓形，半徑為m。由于繞x軸旋轉(zhuǎn)135°時，曲線與繞x旋轉(zhuǎn)45°時的曲線相重合，故圖中只展示繞x軸旋轉(zhuǎn)45°時的曲線。從這些極性圖可以明顯看出，利用視電阻率長短軸的相對大小可以有效地判斷各向異性特征和層理面的傾斜情況。

a．主軸電阻率繞z軸旋轉(zhuǎn)0°、45°、90°、135°；b.主軸電阻率繞x軸旋轉(zhuǎn)0°、45°、90°。圖5 各向異性半空間視電阻率極性圖Fig.5 Polarplots of apparent resistivity for an anisotropic half-space

3.3 含異常體的各向異性圍巖模擬

圖6展示了各向異性圍巖和異常體主軸電阻率繞x和z軸旋轉(zhuǎn)不同角度時的視電阻率極性圖。這里圍巖的主軸電阻率為ρx0=100m，ρy0=400m，ρz0=100m，而異常體的主軸電阻率為ρx1=10m，ρy1=40m，ρz1=10m。對于本文設(shè)計的模型(圖3)，當(dāng)收發(fā)距較小(r=50 m)時，電阻率異常中心位置位于異常體內(nèi)，主要反映異常體的電性特征；而當(dāng)收發(fā)距較大時(r=140 m)，電阻率異常不再關(guān)于異常體中心對稱，主要反映圍巖的電性特征。這為我們識別圍巖和異常體各向異性提供了依據(jù)。

圖6a描述的是異常體和圍巖的主軸電阻率沒有發(fā)生旋轉(zhuǎn)時的視電阻率分布：小收發(fā)距對應(yīng)的視電阻率在x方向被拉伸(長軸)，在y方向被壓縮(短軸)，與圖5中未旋轉(zhuǎn)時的視電阻率分布特征一致；大收發(fā)距對應(yīng)的視電阻率在x方向被拉伸(長軸)，在y方向被壓縮(短軸)，與上述圍巖半空間的視電阻率分布特征一致。圖6b描述的是異常體主軸電阻率繞x軸旋轉(zhuǎn)135°和圍巖主軸電阻率繞z軸旋轉(zhuǎn)45°時對應(yīng)的視電阻率特征分布，與圖6a相比：虛線表示的視電阻率短軸被明顯拉伸，與圖5中主軸電阻率繞x軸旋轉(zhuǎn)135°時的分布特征一致；而實線相對于圖6a沿順時針方向旋轉(zhuǎn)了45°，這與圖5中主軸電阻率繞z軸旋轉(zhuǎn)45°的事實一致。圖6c描述的是異常體主軸電阻率繞x軸旋轉(zhuǎn)90°和圍巖主軸電阻率繞z軸旋轉(zhuǎn)90°時的視電阻率分布，與圖6a相比：由于此時異常體變成水平各向同性，所以虛線短軸被拉伸而成為圓形，與圖5中主軸電阻率繞x軸旋轉(zhuǎn)90°時的分布特征一致；實線電阻率分布相對于圖6a沿順時針方向旋轉(zhuǎn)了90°，這與圖5中主軸電阻率繞z軸旋轉(zhuǎn)90°的情況一致。圖6d描述的是異常體主軸電阻率繞z軸旋轉(zhuǎn)90°和圍巖主軸電阻率繞x軸旋轉(zhuǎn)45°時的視電阻率分布，與圖6a相比：虛線視電阻率沿順時針方向旋轉(zhuǎn)了90°，與圖5中主軸電阻率繞z軸旋轉(zhuǎn)90°的分布特征一致；相比之下，實線視電阻率相對于圖6a短軸被拉伸，下方被拉伸程度大于上方，這是由于異常體存在造成的。圖6e描述的是異常體主軸電阻率繞z軸旋轉(zhuǎn)45°而圍巖主軸電阻率繞x軸旋轉(zhuǎn)90°時的視電阻率分布，與圖6a相比：小收發(fā)距視電阻率繞順時針旋轉(zhuǎn)45°，與圖5中主軸電阻率繞z軸旋轉(zhuǎn)45°一致；而大收發(fā)距視電阻率的短軸被拉伸，呈現(xiàn)圓形分布，與水平各向同性圍巖的視電阻率分布特征一致。圖6f描述的是異常體主軸電阻率繞x軸旋轉(zhuǎn)45°和圍巖主軸電阻率繞x軸旋轉(zhuǎn)135°時的視電阻率分布，與圖6a相比：虛線短軸被拉伸，與圖5中主軸電阻率繞x軸旋轉(zhuǎn)45°時的分布特征一致；而實線短軸被拉伸，由于異常體的存在造成上方被拉伸程度大于下方，與圍巖主軸電阻率繞x軸旋轉(zhuǎn)135°的情況吻合。對比圖6d和圖6f還可以看出：圍巖主軸電阻率繞x軸分別旋轉(zhuǎn)45°和135°時，雖然視電阻率短軸均被拉伸，但由于旋轉(zhuǎn)角度導(dǎo)致各向異性主軸方向不同，拉伸方向相反。這為有效地識別各向異性主軸電阻率旋轉(zhuǎn)特征提供依據(jù)。

旋轉(zhuǎn)角度：a. α0/β0/γ0=0°/0°/0°，α1/β1/γ1=0°/0°/0°；b. α0/β0/γ0=0°/0°/45°，α1/β1/γ1=135°/0°/0°；c. α0/β0/γ0=0°/0°/90°，α1/β1/γ1=90°/0°/0°；d. α0/β0/γ0=45°/0°/0°，α1/β1/γ1=0°/0°/90°；e. α0/β0/γ0=90°/0°/0°，α1/β1/γ1=0°/0°/45°；f. α0/β0/γ0=135°/0°/0°，α1/β1/γ1=45°/0°/0°。圖6 各向異性半空間中存在各向異性異常體的視電阻率極性圖Fig.6 Polarplots of the apparent resistivity for an anisotropic abnormal body embedded in an anisotropic half-space

4 結(jié)論

本文基于非結(jié)構(gòu)網(wǎng)格自適應(yīng)有限單元法，實現(xiàn)了各向異性情況下三維直流電阻率正演模擬，通過與一維層狀介質(zhì)各向異性響應(yīng)對比驗證了該算法的精度和有效性。通過采用環(huán)形掃描測量裝置，對半空間任意各向異性介質(zhì)視電阻率分布特征進(jìn)行分析，并以此為參照對圍巖和異常體各向異性的視電阻率響應(yīng)特征進(jìn)行分析，得出如下結(jié)論：

1)環(huán)形掃面測量方法可以有效識別大地各向異性電性分布特征。較之于傳統(tǒng)的單測線測量的視電阻率，環(huán)形掃面測量結(jié)果可以更好地識別各向異性主軸方向和電阻率等參數(shù)。

2)環(huán)形掃面測量獲得的視電阻率存在各向異性反?，F(xiàn)象，在識別地下各向異性電性分布特征時需要引起注意。

3)利用本文提出的環(huán)形掃面技術(shù)，圍巖和異常體的各向異性特征可以通過改變收發(fā)距分別進(jìn)行識別。

本文對各向異性影響特征的研究和識別方法對于直流電阻率法在淺地表勘查領(lǐng)域的應(yīng)用具有參考意義。

致謝：特別感謝對本文繪圖和算法編寫給予指導(dǎo)和幫助的博士研究生孫思源。

(

)：

[1] 徐世浙，劉斌，阮百堯. 電阻率法中求解異常電位的有限單元法[J]. 地球物理學(xué)報, 1994, 37(2): 511-515.

Xu Shizhe, Liu Bin, Ruan Baiyao.The Finite Element Method for Solving Anomalous Potential for Resistivity Surveys[J]. Chinese Journal of Geophysics, 1994, 37(2): 511-515.

[2]Pain C C, Herwanger J V, Saunders J H, et al. Anisotropic Resistivity Inversion[J]. Inverse Problems, 2003, 19(5): 1081-1111.

[3] 劉斌，李術(shù)才，李樹忱，等. 隧道含水構(gòu)造電阻率法超前探測正演模擬與應(yīng)用[J]. 吉林大學(xué)學(xué)報(地球科學(xué)版), 2012, 42(1):246-253.

Liu Bin, Li Shucai, Li Shuchen, et al. Forward Modeling and Application of Electrical Resistivity Method for Detecting Water-Bearing Structure in Tunnel[J]. Journal of Jilin University (Earth Science Edition), 2012, 42(1): 246-253.

[4]Yin C, Weidelt P. Geoelectrical Fields in a Layered Earth with Arbitrary Anisotropy[J]. Geophysics, 1999, 64(2): 426-434.

[5] 賁放，劉云鶴，黃威，等. 各向異性介質(zhì)中的淺海海洋可控源電磁響應(yīng)特征[J]. 吉林大學(xué)學(xué)報(地球科學(xué)版), 2016, 46(2):581-593.

Ben Fang, Liu Yunhe, Huang Wei, et al. MCSEM Response for Anisotropic Media in Shallow Water[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(2): 581-593.

[6]Yin C. Geoelectrical Inversion for a One-Dimensional Anisotropic Model and Inherent Non-Uniqueness[J]. Geophysical Journal International, 2000, 140(1):11-23.

[7]Habberjam G. Apparent Resistivity, Anisotropy and Strike Measurements[J]. Geophysical Prospecting, 1975, 23(2):211-247.

[8] Li P, Uren N. Analytical Solution for the Point Source Potential in an Anisotropic 3-D Half-Space：Ι: Two- Horizontal-Layer Case[J]. Mathematical and Computer Modelling, 1997, 26(5):9-27.

[9] Yin C, Maurer H M. Electromagnetic Induction in a Layered Earth with Arbitrary Anisotropy[J]. Geophysics, 2001, 66(5): 1405-1416.

[10]Li P, Uren N. The Modelling of Direct Current Electric Potential in an Arbitrarily Anisotropic Half-Space Containing a Conductive 3-D Body[J]. Journal of Applied Geophysics, 1997, 38(1): 57-76.

[11] Wang T, Fang S. 3-D Electromagnetic Anisotropy Modelling Using Finite Difference[J]. Geophysics, 2001, 66(5): 1386-1398.

[12] Wang W, Wu X, Spitzer K. Three-Dimensional DC Anisotropic Resistivity Modelling Using Finite Elements on Unstructured Grids[J]. Geophysical Journal International, 2013, 193(2): 734-746.

[13] Yin C, Zhang P, Cai J. Forward Modelling of Marine DC Resistivity Method for a Layered Anisotropic Earth[J]. Applied Geophysics, 2016, 13(2):279-287.

[14]Si H. A Quality Tetrahedral Mesh Generator and Three-Dimensional Delaunay Triangulator[D]. Berlin: Weierstrass Institute for Applied Analysis and Stochastic, 2006.

[15] 金建銘, 電磁場有限元法[M]. 西安: 西安電子科技大學(xué)出版社, 2014.

Jin Jianming. The Finite Element Method in Electromagnetics[M]. Xi’an: Xi’an University Press, 2014.

[16]Ren Z, Tang J. A Goal-Oriented Adaptive Finite-Element Approach for Multi-Electrode Resistivity System[J]. Geophysical Journal International, 2014, 199(1): 136-145.

[17] 劉國興. 電法勘探原理與方法[M]. 北京：地質(zhì)出版社, 2005.

Liu Guoxing. Principals and Methods of Electrical Prospecting[M]. Beijing: Geological Press, 2005.

[18] Zienkiewicz O C, Zhu J Z. Super-Convergent Patch Recovery and a Posteriori Error Estimates: Part1: The Recovery Technique[J]. Internatioanl Journal for Numerical Methods in Engineering, 1992, 33(7): 1331-1364.

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

Wang Feiyan. 2.5-D DC Resistivity Modeling by the Adaptive Finite-Element Method with Unstructured Triangulation[D]. Changsha: Central South University, 2009.

[20]Wait J R. Current Flow into a Three-Dimensionally Anisotropic Conductor[J]. Radio Science, 1990, 25(5):689-694.