1、二維電導率各向異性介質電磁感應效應的有限元算法李予國 摘要 二維電性各向異性大地電磁場通過有限元方法計算獲得。模型由包含各向異性塊的層狀介質構造組成。每個塊或者層可能是由3×3的電導率張量限定。正演問題可以歸結為兩個耦合的關于垂直-平行于場分量的和。它們由有限元法數值地解出。線性有限元的解系是用預條件共軛梯度法求取的。之后，地表垂直-平行的場分量和，利用樣條插值通過對和進行數值微分得到。 二維有限元算法通過對比二維有限差分的解證明是正確的。三個類型的模型的大地電磁響應用來證明各向異性的影響：水平，垂直，傾斜各向異性。第四個模型用來模擬在剪切和俯沖帶各向異性的影響。這些模型的響應模擬了

2、電阻率曲線在長周期和具有明顯的主對角線元素的張量主抗存在室的分離，正如以前觀測到的。1、簡介 近些年來，越來越多的關注投入到了對電性各向異性的電磁感應方面上的研究，尤其是為了試圖完全理解對更長周期的大地電磁的觀測數據。加拿大地區的大地電磁測量顯示了電性的各向異性在下地殼和上地幔（Kellett等，1992；Mareschal等1995）。大的來自圍繞德國深鉆點（KTB）各向異性大地電磁曲線解釋為上地殼高度的電各向異性（Eisel 和 Haak 1999）。Rasmussen (1988)使用一個深地殼層內的各向異性模型來解釋瑞典南部的大地電磁橫斷面數據。層狀構造的各向異性效應首先被OBrien

3、 和Morrison (1967)研究。Reddy 和 Rankin (1975)首先開始研究二維各向異性模型，他們只考慮水平各向異性的工作。最近，Osella和Martinelli (1993)計算了一些圓滑不規則邊界有一定的主軸旋轉角的大地電磁響應。Schmucker (1994)提出一個計算不均勻薄的覆蓋體在層狀半空間的電磁感應效應，其中可以含有一個或者兩個各向異性的電導層。用有限元方法，Pek 和 Verner (1997) 還有 Weidelt (1996) 模擬了廣義二維、三維各向異性構造，分別具有任意的主軸旋轉角。 本文中，二維模型被重新研究，但是使用有限元（FE）的方法。首先，

4、我們詳細描述有限元程序的數值近似。之后，我們展示有限元程序模擬不同的簡單測試模型的響應。我們的結果和Pek& Verner (1997)的有限差分程序進行了對比。最后，我們計算三種類型的各向異性介質的大地電磁響應：水平，垂直，以及傾斜各向異性。我們最后以模擬一個地質剪切帶或者板塊俯沖帶的地質背景下的模型結束。本文這樣的構造是為了展示一個有限元策略去實現二維各向異性介質的電磁感應。盡管有限差分對于這類問題是有效的，有限元解仍然是需要的，因為有限元和有限差分都有自己的特殊的優勢還可以彼此互相核對。另外，有限元方法可以模擬非矩形線性構造的真實地球構造。2、邊界值問題 考慮到圖1中的二維模型。

5、陰影異常區嵌入到一個簡單的層狀構造中，這個構造包括n層而且底部的層延伸到無窮遠。為了簡單，異常區域沒有顯示內邊界作為一個定向獨立的電導率，它至少和周邊的第j層的電導率張量不同。然后我們的程序允許更一般的模型，異常區域可以細分成不同的不均勻塊而且可以接觸不同的層。 異常區域延伸到y正向到無窮遠，例如模型可以出現到一個新的正常構造。在這種狀態下算法排除源的感應。模型在走向方向x方向是不變的。感應電磁場在x方向也是不變的，盡管如此場的矢量也有三個分量對于極化的磁源平行或者垂直走向，由于各向異性。因此，通常區分TE和TM模式的各向同行構造就變得不合理了。由此就有結果，當主磁場矢量垂直軸向，電場將有主場

6、方向的分量，導致電荷產生于界面（除了傾斜各向異性的例子）。 假設一個時諧變量，在準靜態近似的電磁場的情況下 (1) 其中是自由空間的磁導率，且 為電性電導率張量。它是對稱的，當旋轉到主軸方向（x',y',z'）上去，可以表示為 圖1.本文中考慮的二維各向異性模型在特別的二維例子中，方程（1）消退為常數電導率的均勻塊 很明顯如果垂直-平行分量和得到了，剩下的分量和，和可由和的空間導數得到。方程（2）-（7）可以得到兩個關于分量和的二階偏微分方程： 從方程（8）和（9）很明顯看出各向異性和獨立的垂直-平行分量和通過一階偏微分方程耦合在一起。因此，在這種情況下不存在獨立的TE

7、、TM模式對于廣義各向異性情況。因此，這些方程必須同時利用和求解。 之后敘述了不同的空間形式的各向異性的模擬結果。如圖1一個均勻的異常區域具有各向同性的正常結構的圍巖，一個或者兩個平行于參考主軸（x，y，z）坐標系的主軸被考慮了進去。2.1 水平各向異性 當，主軸Z'是垂直的，另外兩個主軸x'和y'在水平面（x,y）具有相對于x軸的垂直角。方程（8）和（9）就消退為 感應方程變得比方程（8）和（9）更加簡單，但是垂直-平行分量和仍然耦合在一起通過關于z的一階偏導。所以方程（10）和（11）必須聯立在一起去為了得到和。2.2 傾斜各向異性 當，電導率張量主軸x'是

8、水平的在走向方向，另外兩個主軸y'和z'在垂直面（y,z）具有相對于y軸的傾斜角。現在方程（8）和（9）就解耦為兩個獨立的模式方程（12）對應于E極化模式可以通過二維各向同行程序解得-只需要電導率替換為作為標量電導率。然而，TM模式的表達式子仍然比較復雜。2.3 垂直各向異性 如果，電導率張量的三個主軸真好和坐標系（x，y，z）是一致的。如果再有，這種結果關于z軸軸對稱，水平方向的電導率完全不同于垂直方向的電導率。這樣（8）和（9）仍然是耦合在一起的，我們有 盡管B極化模式感應方程比方程（9）更加簡單，它不能通過二維各向同性程序解得，因為電導率在水平方向和垂直方向是不同的。如果

9、兩個方向的電導率相同，方程（15）消退為各向同性的例子。之后的邊界條件應用：模型的外邊界條件，設定為狄里克雷邊界條件，由對應的一維層狀模型的左右邊界條件。在模型的頂部和底部由左右邊界的值進行線性插值得到（Pek 和 Verner 1997）。在內邊界上，切向電場和磁場的分量，和必須連續。根據圖1，和為 重構方程式子之間的關系至關重要。設n為不均勻區域的外側的單位矢量，且有 和來自方程（8）和（9），和為y軸和z軸的單位向量相應的。這樣，實質性的偏導數后，表達式為和為發現的切向分量。方程（17）是自動滿足的，對自后的偏導是很重要的。結合這些表達式基本的積分-微分方程（23）和（24）可以得到濃縮

10、的形式，這些形式可以直接正演數值模擬對待。更多細節看李（2000）。3 有限元方法 由方程（8）和（9）組成的數值近視問題，基于有限元方法。垂直-水平分量和跟電導構造相切，所以處處連續。于是，在有限元中假設所有分量穿過單元都是滿足連續。這種近似建立在模型域，整天包圍二維不均勻體，區域延伸的足夠遠使得異常場在邊界處很小。為了避免絕緣的空氣層在地面上產生奇異性，我們假設空氣層具有很小的電導率，但是非零，實際上小于在我們的計算中。數值計算結果顯示電導率到。 加權余量方法用來獲得從偏微分方程（8）和（9）得到的積分方程。方程(8)乘以一個關于電場的任意小的的變量，在模型域里的積分： 上述方程中，第一項

11、包含二階偏導數，可以用格林函數簡化。 其中指示模型域的邊界。之后方程（19）可以寫成等價形式 類似地，方程（9）乘以一個磁場分量任意小的變量,之后用高斯公式改變一下 加上方程（16）這樣得到積分方程 這里我們仍然用這個公式模型域可被分割為矩形或者是三角形單元。矩形單元已被李（2000）年描述過。以下的部分，展示三角形單元。公式（20）和（21）分解到每個單元帶有指示標e=1,2,ne其中表示為三角形單元e的邊界。在形成方程（23）的時候，我們利用方程（18）時替換被積函數為.用一種簡單的方法方程（24）由方程（17）獲得。 現在邊界條件必須被利用。在內部邊界切向電場和磁場（）是被要求的。因為每

12、個邊界在集成的過程中被遍歷兩次且方向相反，則在單元邊界的線積分的總和為零。在外邊界，因為設置的狄利克雷邊界條件，電場和磁場的變量和等于0.因此鮮雞粉也等于0.這樣以來，方程（23）和（24）最終消退為根據有限元的線性近似方法，我們假設在每個三角形單元內，和都為關于y,z的線性函數：其中和為第i個頂點的電場和磁場對應的坐標為（），i=1,2,3為三角形（圖2），為線性基函數。它們由一下定義： 其中 方程（25）和（26）在一個單元內的面積分用方程（27）-（32）估算。所有單元的積分都可以組裝成兩個線性方程系。詳細內容見附錄A。聯立這些方程系，我們可以最終寫下方程以矩陣的形式： 其中 和是階方陣

13、（是結點的數目），和不對稱的方陣，滿足因此矩陣K是一個2的對稱方陣。它是稀疏的復元素。U是2列向量，它是包含電場和磁場的未知數。方程（33）代入外邊界條件，在內部的點上和可以用共軛梯度方法解出來。圖2三角形單元切面 共軛梯度法最早由Hestens 和Stiefel (1952)提出，只適合求實數對稱的正定的矩陣。Jacobs (1981)提出一種B共軛梯度法使得共軛梯度法應用至復矩陣非正定系。廣義的共軛梯度法通過對系數矩陣k的預處理得到改善。在我們的例子中簡單的雅克比擴展是很好的預處理。本文中我們估算的模型，網格包含2500多個點，方程的數目大約5000.我們總是用0向量U=0作為初始值，我們

14、總是采用趨于一致的解。 在共軛梯度方法中，僅稀疏對稱的系數矩陣對角線上和下的元素是需要的，它是儲存到一維數組中根據Schwarz (1991)的方法。利用這種共軛梯度法求解，計算內存大大降低。3 有限元方法 求解線性方程組（23），我們獲得和在每個節點上的值。另外的場分量和很容易從方程（3）和（4）中獲得，而分量和是根據（6）、（7）中一些額外的計算求導獲得的。總結為 在我們的算法中，需要的導數是用樣條插值數值計算得到的。 阻抗張量元素可以用兩個正交的線性源極化模式的電場和磁場計算得來。例如主磁場方向在（極化模式一）。或者在y方向上（極化模式2）.主場為z<0處源的感應場和正常構造在處的

15、誘導場的和。 我們得到阻抗張量元素 電阻率和阻抗相位為 5 數值測試 為了測試電性各向異性有限元的適應性，我們的計算結果對于兩種模式和已經獲得的有限差分方法進行了比較。圖3所示第一個測試模型首先由Reddy &和Rankin (1975)提出。我們計算這個模型在10s時候的響應，且和Pek 和Verner (1997)的計算結果對比。結果展示在圖4中，有限元和有限差分吻合很好。圖3水平各向異性巖脈模型Reddy &和Rankin (1975) 圖4圖三中的視電阻率（上部）和相位（下部）曲線。方塊為Pek 和 Verner (1997)的有限差分結果，實線，為本文描素的有限元結果

16、。圖5 有限元和有限差分測試模型的對比，模型包括一個出露地表的巖石塊，覆蓋在一個各向異性層上面。圖6 圖5中的模型的視電阻率（上部）和相位（下部）圖。方塊為Pek 和 Verner (1997)的有限差分結果，實線，為本文描素的有限元結果。 圖5展示第二個模型有Pek 和 Verner (1997)年提出。一個水平各向異性層位于一個二維廚樓地表的各向異性場下。所涉及到的兩個構造的走向彼此平行，但是不平行于二維模型的走向。模型可以用來展示MT數據的擾動，可以用一個復雜的各向異性例子計算得到。在圖6中，有限元計算結果對于30s和Pek & Verner (1997)的對比。電阻率的相對誤差

17、低于0.5%，阻抗的差異不超過1°。6 各向異性效應 已經注意橫向非均勻各向同性地球模型不足以構成大地電磁的阻抗張量，特別是各向 異性阻抗保持或多或少的不變在延伸區域。這在長周期已經很好的研究，最好在一天之內的變化的周期，曲線最大值和最小值能差兩個數量級或更多。的最小值在10（Schmucker 1998）。在結晶巖的表面，類似的觀測已經到短周期內，我們的模型已應用到。響應的在x'主軸電阻率我們選擇100-500，x'主軸電阻,10,z'主軸電阻率100-500. 在這部分里，一個簡單的巖脈模型，在一個各向同性的半空間內，用來展示水平，垂直和傾斜各向異性的大地

18、電磁響應。所有的計算都在t=10s，除非有另外的敘述。對于均勻半空間10s的屈服深度為50km，5km對于10在y的主軸方向。值得一提的是維度的板狀體，低邊界在9km比半空間的屈服深度小很多。6.1 傾斜各向異性圖7 一個具有傾斜異常的巖脈在各向同性的均勻半空間。巖脈的主軸電阻為，且具有變化的傾向角。圖8 圖7中10s時不同的變化角視電阻率圖圖9 顯示圖7中的電場在yz平面內的分布。產生的主場為，計算的周期為10s。 圖7展示主軸面向切面的巖脈。各向異性巖脈的電阻率分別為。圖8顯示視電阻率對于變化的角在10s。這張圖表明： （1）視電阻率獨立于傾角，大地電磁場僅單向取決于。在這種情況下，且不受

19、各向異性的影響。 （2）電阻率受的影響。曲線不再對稱關于模型的中心。模型的最小值偏離中心，偏向的一側取決于傾角的一側。偏離會變大隨著水平或者垂直方向的角度偏差變大。我們解釋這種現象通過看電場的分布。圖9展示電場的分布在yz平面對于三個傾角（0°,30°,60°,90°）在極化模式1，產生主場的方向是。箭頭的起點是有限元點。箭頭長度和方向代表節點處場的大小和方向。圖很明顯地顯示場的大小和方向變化隨著各向異性巖脈傾角變化。（3）=0°時候，對應于一個y方向10，x和z方向500的模型響應。類似的=90°，對于500在x，y方向，z方向10

20、.視電阻率和相位對應于變化的角（0°,30°,60°,90°）在中心c（y=0）在圖（a）和（b）分別展示。視電阻率在非常低的周期接近均勻本空間的真實電阻率，相位接近45°。然而，隨著周期增長，和分離很明顯，和同樣，偏離大小取決于傾角。然而和因不同而不同，和曲線在長周期也適合均勻半空間的電阻率，相位和近似45°，但仍然受各向異性巖脈的影響。關于c點對稱的a點和b點的和在10（c）和10（d）上分別顯示。從這些圖中可以看出來，視電阻率和相位曲線在a點，和不同于b點。然而當=0或90°時候，卻是相同的。圖10 圖7中A、B、C地

21、表點的視電阻率和相位曲線。6.2 水平各向異性 圖11 一個2-D水平各向異性巖脈在均勻各向同性圍巖中。巖脈的主軸電阻為，且具有變化的各向異性走向角。圖12 大地電磁主軸具有變化的走向角。叉表示值大地電磁張量非對角元素和，棒代表主軸元素的值和，參考Swift-rotated坐標系（x''，y''）。圖13電場實部Re（Ex）和Re（Ey），在模型11的表面。產生的主場為，計算的周期為10s。 圖11顯示了主軸方向，它們現在水平方向傾斜角。不均勻的主軸電阻率分別。圖12顯示大地電場阻抗用西門子的表述方法沿著一個剖面表面，對于變化的各向異性角在10s。從圖12我們可

22、以得到： （1）各向異性塊上。最大和最小軸顯示最大和最小的電導率。 （2）張量點到率主對角線元素和存在（除了30°和90°），且隨著角增加而增加。 （3）遠離不均勻體，異常場消退，非對角線元素變的相等和，對角元素消失：。當時。當從各向異性場通過各向同性場的時候，阻抗會在過渡區有一定的減少，這里的方向可能會嚴重扭曲，特別是涉及到各向異性體和淺部異常體。圖13顯示電場實部Re（Ex）和Re（Ey），在剖面的表面具有不同的角（=0°，30°，45°，60°，90°）在第二種極化模式，產生的磁場。可以從圖中看出隨著各向異性角的改變。

23、場的大小和方向都在改變，y方向電場，Re（Ey）存在=0°或=90°；當在兩個角之間時Re（Ey）變大。6.3 傾斜巖脈模型 圖14 圖7中的一個二維模型=90°，巖脈旋轉了45°模擬方向獨立的剪切和俯沖帶，具有不同的平行和垂直電阻率。 圖14顯示了一個圖7中的二維模型=90°，但是巖脈旋轉了45°，模擬方向獨立的剪切和俯沖帶。由于剪切作用，巖脈在剪切面和垂向有不同的電導率。如果沒有不同存在，假設整個巖脈作為一個整體比圍巖的導電性好，一個關于中心不對稱的最小值在走向方向被觀測到。如果平行電導率10倍阻性圍巖。對稱性得到了保證單在巖脈上

24、最小值不可見。如果再剪切方向更加導電，最小值解更加顯著，正如期待的（圖15a），關于分歧。圖15（b），1s一下開始，在周期更長時編程一個靜的偏移。圖15視電阻率在巖脈的剖面上對于10s圖14中的模型（a）,測量點位于y=0（b）。電阻率垂直于巖脈邊緣（各向同性例子），或者十倍小于。視電阻率曲線對于在10-1000之間變化，曲線緊挨著。7 結論 我們展示了一個二維大地電磁數值迷你程序。其主要特點如下：（1）每個各向異性塊的電導率張量都由一個3乘3的矩陣構成，這樣保證任意的各向異性構造都能被考慮，包括空間中的水平、垂直、傾斜的各向異性。（2）數值近似問題由有限元方法得到，這種方法很適合具有坡度的

25、邊界模型和地形，如圖14中所舉的例子。（3）假設一個非常小，但正的空氣介質，空氣層可以和導電模型成為一個整體。兩種模式的方程在整個模型區域均勻近似，這樣就簡單一些相對于區分空氣和大地或者僅大地變量如Pek 和 Verner(1997)，它們把空氣的電導率設為0.線性有限元方程利用了快速迭代的技術，增加的空氣層中的變量Hx不是過多的問題。數值試驗顯示場分量是穩定的，實際上不熟寬的空氣電導率的影響限制，特殊的從10-12開始，這對偃師的電導率已經足夠低，已經低于大地上面的空氣電導率。 （4）方程（33）是利用改變的共軛梯度法求的，是一個復的對稱矩陣。我們發現簡單的雅克比擴展是一種有效的預處理方法。

26、 （5）垂直走向場Ey和Hy用通過樣條插樣數值得到。致謝 作者對U. Schmucker博士和J. Pek博士對本文的討論和改進表示感謝。J. Pek博士也提供了MT一維各向異性構造的正演程序和兩個二維的有限差分的結果。同時也感謝Martyn Unsworth, Chester Weiss和另外一位匿名的編輯的建設性意見。作者感謝德國交流服務提供的博士獎學金。本文由格延根大學地球物理許愿管理。參考文獻 Eisel, M. & Haak, V., 1999. Macro-anisotropy of the electrical conductivity of the crust: a m

27、agneto-telluric study from the German Continental Deep Drilling site (KTB), Geophys. J. Int., 136, 109122. Hestens, M.R. & Stiefel, E., 1952. Method of conjugate gradients for solving linear systems, J. Res. Nat. Bur. Stand., 49, 409436. Jacobs, D.A.H., 1981. The exploitation of sparsity by iter

28、ative methods,in Sparse Matrics and their Uses, pp. 191222, ed. Duff, I.S.,Springer, Berlin. Kellett, R.L., Mareschal, M. & Kurtz, R.D., 1992. A model of lower crustal electrical anisotropiy for the Pontiac Subprovince of the Canadian Shield, Geophys. J. Int., 111, 141150. Li, Y., 2000. Finite e

29、lement modeling of electromagnetic elds in two- and three-dimensional anisotropic conductivity structures, PhD thesis, University of Go ¨ ttingen. Mareschal, M., Kellett, R.L., Kurtz, R.D., Ludden, J.N., Ji, S.,& Bailey, R.C., 1995. Archaean cratonic roots, mantle shear zones and deep elect

30、rical anisotropy, Nature, 375, 6527, 134137. OBrien, D.P. & Morrison, H.F., 1967. Electromagnetic elds in an Nlayer anisotropic half space, Geophysics, 32, 668677. Osella, A.M. & Martinelli, P., 1993. Magnetotelluric response of anisotropic 2-D structures, Geophys. J. Int., 115, 819828. Pek,

31、 J. & Verner, T., 1997. Finite-difference modelling of magneto-telluric elds in two-dimensional anisotropic media, Geophys. J. Int.,128, 505521. Rasmussen, T.M., 1988. Magnetotellurics in southwestern Sweden;evidence for electrical anisotropy in the lower crust?, J. geophys. Res.,93, 78977907. Reddy, I.K. & Rankin, D., 1975. Magnetotelluric response of laterally inhomogeneous and anisotropic media, Geophysics, 40, 10351045. Schwarz, H.R., 199