基于環形掃面測量的三維直流電阻率法任意各向異性模型響應特征

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

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

0 引言

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

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

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

1 正演理論

1.1 電導率張量

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

(1)

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

σ=RTσcR。

(2)

式中：σc為主電導率張量，

(3)

(4)

式中，α、β、γ分別是電導率主軸繞x、y、z軸的旋轉角[6]。

1.2 基于二次場的邊值問題

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

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

(5)

(6)

σ·up·n-

(7)

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

(8)

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

?Ω(·(σ·us)+

(9)

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

KUs=b。

(10)

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

1.3 基于梯度恢復的誤差估計和局部網格自適應加密策略

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

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

i=1, 2,,n。

(11)

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

du=Ga。

(12)

式中：

(13)

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

在自適應網格加密策略中，我們首先設定每個目標區域的最大誤差閾值:

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

(14)

式中，k表示目標區域數量。進而，可以求出每個單元一次迭代時的最大體積：

(15)

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

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

2 精度驗證

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

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

通過分析圖2可以得出如下結論：

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

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

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

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

3 數值實驗

3.1 模型建立

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

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

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

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

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

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

3.2 各向異性半空間模擬

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

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

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

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

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

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

旋轉角度：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 結論

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

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

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

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

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

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

(

)：

[1] 徐世浙，劉斌，阮百堯. 電阻率法中求解異常電位的有限單元法[J]. 地球物理學報, 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] 劉斌，李術才，李樹忱，等. 隧道含水構造電阻率法超前探測正演模擬與應用[J]. 吉林大學學報(地球科學版), 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] 賁放，劉云鶴，黃威，等. 各向異性介質中的淺海海洋可控源電磁響應特征[J]. 吉林大學學報(地球科學版), 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]. 西安: 西安電子科技大學出版社, 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]. 北京：地質出版社, 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] 王飛燕. 基于非結構化網格的2.5-D直流電阻率法自適應有限元數值模擬[D]. 長沙：中南大學, 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.