等效場源法的CSAMT三維無限元正演模擬

張林成 胡宏伶 湯井田 肖衛初 肖 曉 原 源

( ①湖南城市學院信息與電子工程學院，湖南益陽 413002；②中南大學地球科學與信息物理學院湖南長沙 410083；③計算與隨機數學教育部重點實驗室(湖南師范大學數學與統計學院)，湖南長沙 410081；④核資源與環境國家重點實驗室(東華理工大學)，江西南昌 330013 )

0 引言

隨著深部資源勘探的不斷深入，要求可探測深度越來越大[1]。可控源音頻大地電磁測深法 (CSAMT，Controlled Source Audio Electromagne-tic Method)[2]以勘探深度大、抗干擾能力強、工作效率高等特點成為研究熱點，是中深部資源勘探的重要手段之一。該方法在隱伏礦體勘察、復雜地區頁巖氣勘探、地熱資源調查及海洋油氣勘探等領域都取得了良好效果[3-4]，成為發展快且前景可期的地球物理勘探方法之一。

CSAMT電磁場正演的主流方法主要有邊界單元法、有限差分法、積分方程法和有限單元法四種，其中有限單元法以其理論系統化、適應性強、計算精度高、弱解可微等優點得到了更多的重視和應用。Coggon[5]于1971年提出了大地電磁問題的有限元模擬方法，自此有限元計算在電磁領域得到了極大的發展。CSAMT有限元數值模擬中，邊值問題包括控制方程和邊界條件兩個方面。對于控制方程，由于場源存在奇異性，場源的處理方式是關鍵，常用的方法有二次場法和總場法，其中二次場法是主流。二次場法將場分解為背景場和異常場，背景場利用均勻半空間或層狀模型解析解可直接計算，二次場則通過有限元法求取[6-12]。總場法直接從總場著手，采用近似法模擬奇異性場源特征(例如偽delta函數法)，然后通過有限元求解場值[13-22]。在CSAMT三維有限元正演模擬中，無論采用總場法還是二次場法，場值的求解精度都是正演模擬是否成功的標志，因此開展不同場源方案下場值求解精度的研究十分必要。

一般來說，對于邊界條件的處理基本上采用傳統截斷邊界方法[6-22]，即在一個相對較大的區域內，將無窮邊界問題近似為有限區域，這往往會造成有限元計算區域太大、節點數太多、存儲量過大和計算耗時過長等問題。較合理的邊界處理策略應該是將計算區域限制在較小的目標區域內，以減少單元數和節點數，實現節約計算資源及提高計算速度的目的。有限元-無限元耦合法是針對常規截斷邊界引起的大尺度模型問題提出的[23]，具有離散區域小、求解速度快和計算精度高等優點。其核心思路是用“無限單元”取代傳統截斷邊界，通過坐標映射將“無限單元”沿某一方向無限擴展，實現電磁場在無限單元中快速衰減為零，最后通過有限元-無限元耦合求解方程組，從而實現CSAMT的三維正演計算。

在地球物理領域，無限元思想主要應用于地震、測井和直流電領域。Blome等[24]應用Astley無限元處理直流電法邊界問題，縮小了計算范圍，減少了節點數, 且避免了采用Dirichlet或混合邊界條件，但是其采用的無限元形函數較復雜，且由于引入了Astley無限元中特有的額外權重因子，使整體系數矩陣失去了對稱性。Nath等[25-26]將無限元分別應用于重磁法的數值計算。在電磁法領域，公勁喆[27]將無限元應用到直流電阻率法的三維正演模擬。湯井田等[28]將基于有限元-無限元的三維直流電成功應用于泥河鐵礦勘探實例；張林成等[10]、湯井田等[22]采用無限元代替截斷邊界條件，成功地將有限元-無限元耦合法應用于CSAMT法，實現了小區域內CSAMT數據的三維快速高精度正演，但其采用的是基于電場二次場的方案。

針對以上研究現狀，本文從基于總場方案的邊值問題出發，首先采用等效場源法，假設場源周圍很小區域內電磁場的分布是均勻的，用區域內網格節點上電磁場的近場值表示場源，并將節點場值作為第一類邊界條件帶入場方程，實現場源的近似模擬；然后，用“無限單元”代替傳統截斷邊界條件，采用Garlerkin有限元法和Pardiso并行直接求解器實現了三維CSAMT快速正演計算；最后，通過均勻半空間模型解析解驗證方法的正確性，并開展了場源等效模擬最佳范圍的研究。

1 邊值問題

1.1 基本方程

對于角頻率為ω、時間因子為e-iω的水平電偶極子，在各向同性導電介質中產生的電場E滿足雙旋度方程[29]

(1)

式中：σ為電阻率；Js為電偶極子電流密度矢量；μ0表示真空磁導率，取空氣和地下介質的磁導率均等于μ0；ε0是真空介電常數。

1.2 場源的計算

由于水平電偶極源存在奇異性，如何精確加載場源成為基于電場總場數值模擬的關鍵問題。本文采用等效場源法，將其周圍很小的區域看成均勻的，用此區域內離散網格節點上電磁場的近場值表示近似等效場源，并將源附近節點上的場值第一類邊界條件代入方程。最后通過有限元—無限元耦合求解方程組，從而實現三維正演計算。

設在無窮大的空間中，存在一無限大分界平面S，把無窮空間分為上下兩部分，如圖1所示。在分界面S上h高度處有一水平電偶極子AB，其偶極距為dL。建立直角坐標系和柱坐標系，選取共同原點，位于偶極子中心，偶極距與x軸正方向一致，z軸垂直向下。設上半空間和下半空間介質的電導率分別為σ1和σ2，導磁率為μ1和μ2，介電常數為ε1和ε2。在勘探地球物理中，通常下半空間代表大地，上半空間代表地面之上的空間。設偶極源中的電流為正弦交變電流I=I0e-iωt，其中I0為電流振幅。

在柱坐標系中，上、下半空間中任意一點的電場分量為

(2)

(3)

根據電法理論，在電性分界面上，電場的法向分量是不連續的，而節點有限元要求電磁場法在法向上是連續的。因此，求解得到的有限元解往往不準確，需要對解進行校正。

在有源區，節點有限元的電場解不滿足

·(εE)=-·[JS/(iω)]

在無源區域不滿足

·(εE)=0

因此，需在式(1)中加入散度校正條件[25]

(4)

圖1 均勻半空間參數及坐標示意圖

1.3 Galerkin有限元法

基于加權余值有限元法[30]將場源作為第一類邊界條件，建立式(4)的殘差公式

r=××E+k2E+·E

(5)

取任意測試函數V，在區域Ω上令

?Ωr·VdΩ=0

有

(6)

設T為區域Ω的外邊界面。用無限元處理無窮邊界問題，在無限單元內電磁場衰減為零，即

E=0

利用第一矢量Green定理[31]可將式(6)簡化為

(7)

利用有限元對內部計算區域進行離散。假設內部計算區域內有n個節點，對于第j個測試函數Vj有

(8)

2 有限元—無限元耦合算法

有限元—無限元耦合算法即將整個求解區域分為有限元區域和無限元區域，用無限元區域代替傳統的外邊界，在這兩個區域分別采取有限元分析和無限元分析，通過總體剛度矩陣組裝將兩者結合，從而進行數值求解。

圖2為有限元和無限元計算區域劃分示意圖。圖中灰黑色區域為有限元區域，即目標區域，包含場源、目標體及測點等；白色區域為無限元區域，即從有限元區域邊界至無窮遠，是邊界計算區域。無限元分析就是在該區域的某一方向上采用無限元映射和形函數，將整體坐標映射到局部坐標上，其原理與有限元分析相同。

圖2 有限元和無限元計算區域劃分示意圖

圖中P為坐標原點，節點1～8為無限元基本要素，節點1～4為有限元邊界上某單元上的節點，P點到節點5～8的距離分別為P點到節點1～4距離的兩倍，圖4同

2.1 有限元分析

有限元單元分析時，采用矩形六面體進行區域離散，單元節點編號及坐標如圖3所示。

圖3 有限元映射示意圖 (a)子單元； (b)母單元

圖3中，子、母單元坐標對應關系為

(9)

矩形六面體形函數表達式為

(10)

式中(ξq，ηq，ζq)(q=1，2，…，8)表示子單元節點q在母單元中的坐標。

2.2 無限元分析

無限元單元分析時，采用三維八節點Astley型無限元。圖4是三維無限元映射示意圖。無限元分析就是通過無限元映射將無窮坐標映射到圖4b中的局部坐標上。圖4b中無限單元最外圍的四個節點9～12代表無窮遠，其場值為零。

圖4 無限元映射示意圖 (a)子單元； (b)母單元

圖4b中ξ方向表示無限元映射方向，在ζ-η平面內，無限單元與有限單元采用相同的映射形式和形函數。無限元坐標映射為

(11)

再結合ξ方向的坐標映射關系式，可得到8個結點無限單元坐標映射函數

(12)

無限元形函數M的表達式為

(13)

上式即Astley映射無限單元理論[30]中采用的形函數中的系數(1-ξ)/2與二階Lagrange插值多項式的乘積。

有限元與無限元單元分析基本一致，都為八節點單元，因此在數值模擬中，可將無限元和有限元完美結合起來，保證剛度矩陣的對稱性，從而使得求解更簡單、更方便。

3 方程組求解

三維CSAMT有限元-無限元耦合正演模擬最終形成大型、稀疏、對稱的復系數線性方程組

Ax=b

(14)

式中：矩陣A為混合有限元-無限元總體剛度矩陣，為3×Nx×Ny×Nz階方陣，其中Nx、Ny、Nz分別為x、y、z三個方向上的節點數，矩陣A中每行最多有81個非零元素；x為各節點待求解的電場值；b為右端項。對于大型稀疏線性方程組，本文采用性能良好、高度并行化的開源求解器Pardiso, Pardiso直接求解，采用LU分解，這種方法特別適合于多源CSAMT的情況。

4 精度驗證

為驗證本文基于等效場源的CSAMT三維無限元數值模擬程序的正確性，對本文方法計算結果與均勻半空間地表水平電偶極子產生的電磁場解析解進行對比、分析。以下數值計算平臺均為Intel(R) Xeon(R) CPU 3.10G，256GB RAM，16CPUs。

4.1 模型計算

假設一個均勻半空間模型：地下介質電阻率為100Ω·m，空氣中電阻率為1×109Ω·m，水平電偶極子沿x軸布設，長度為1m，供電電流為1A，水平電偶極子中心點位于坐標原點，計算頻率為256Hz。對模擬區域進行網格剖分，x方向區域為-5000～5000m，網格距為100m,共計101個網格節點，左右對稱；y方向剖分方案與x方向一致；z方向(空氣中z坐標為負)區域為-2000～2000m，網格距不等，場源附近較小，共35個網格節點，地面以上部分剖分為9層。

圖5為均勻半空間模型下基于等效場源的CSAMT三維有限元-無限元正演模擬數值解與解析解對比。場源加載范圍為1.5倍趨膚深度, 所示結果為z=0、y=1500m、x為-1000～1000m剖面，頻率為256Hz時，電場Ex和Ey振幅數值解與解析解及相對誤差。

從圖5a可以看出，基于等效場源的CSAMT三維無限元數值計算的電場Ex和Ey數值解與解析解場值大小非常相近，且曲線變化形態一致。圖5b中，Ex數值解相對解析解最大誤差為1.5%，平均相對誤差小于0.8%；Ey數值解相對解析解誤差較Ex稍大，但均小于2.8%，平均相對誤差約為1.5%。因此，基于等效場源的CSAMT三維無限元數值模擬程序的計算結果是正確可靠的，其計算精度較高。

4.2 相同區域無限元法與傳統有限元法結果對比

針對上述均勻半空間模型，對相同區域進行相同的網格剖分，利用傳統有限法進行模型正演，并對比分析本文方法與傳統有限元法的計算結果。圖6為無限元法與傳統有限元法電場分量Ex和Ey振幅數值解與解析解的相對誤差。

圖5 均勻半空間模型基于等效場源的CSAMT三維有限元—無限元正演模擬數值解和解析解(a)及相對誤差(b)

由圖6可見，在相同的小區域內，即x為-1000～1000m范圍內，傳統有限元法計算所得電場分量Ex和Ey的平均相對誤差約為5%，遠大于無限元法(小于3%)，即無限元單元法在較小的計算區域內的計算結果精度高于傳統有限元法。因此，本文基于等效場源的CSAMT三維無限元法相對傳統大區域有限元法，能在保證精度的前提下，縮小計算區域，提高計算速度。

4.3 無限元法與傳統大區域有限元法計算效率對比

關于無限元法與傳統大區域有限元法計算效率對比，文獻[10]有詳細分析，此處不再贅述。從表1可以看出，相對于基于總場法的傳統大區域有限元法，本文提出的基于等效場源的CSAMT三維無限元計算性能更高，明顯具有離散區域小、計算節點少、求解速度快等優勢。

圖6 無限元法與傳統有限元法Ex(a)和Ey(b)數值解相對誤差

表1 基于總場法的傳統大區域有限元法與小區域無限元法計算數據對比[10]

4.4 無限元法多頻點計算結果對比分析

針對4.1節均勻半空間模型，分析不同頻率下的計算結果的精度。圖7是頻率分別為64、256、1024Hz時利用本文方法計算的z=0、y=1500m剖面的Ey值及Ey數值解相對誤差。可以看出，在相同剖分方案下，1024、256、64Hz下的電場分量Ex和Ey的數值解平均相對誤差都小于2%，精度較高。還可以看出，隨著計算頻率的降低，本文方法計算的Ex和Ey數值解的計算相對誤差逐漸減小，這是由于頻率降低，趨膚深度增大，均勻半空間模型剖分網格相對加密。因此在進行正演時，為提高精度，可根據頻率的不同，進行不同的網格剖分。

圖7 不同頻率下Ex(a)和Ey(b)數值解相對誤差

5 場源等效模擬最佳范圍研究

等效場源法是將場源周圍很小的區域看成均勻的，用此區域內網格節點上電磁場的近場值表示源，將源附近節點上的場值第一類邊界條件帶入方程。對于水平電偶極源，由于電磁場在近區、過渡區和遠區的衰減規律不一樣，利用節點電場值代替場源時，一方面僅限于近區，另一方面該區域中的網格應足夠多，使得場值能充分模擬場源的特征。因此，如何精確確定源的加載范圍是關鍵。本節以參照CSAMT的近區范圍，通過數值模擬，分析場源不同加載范圍的結果。

圖8所示是場源不同加載范圍對比試驗結果。采用均勻半空間模型，電阻率為100Ω·m，水平電偶極子沿x軸布設，其中心位于坐標原點，長度為1m，供電電流為1A，計算頻率為256Hz。網格剖分與4.1節模型剖分方案一致。計算可知，256Hz對應的趨膚深度δ約為160m，分別對比觀測點到原點的距離r=0.5δ、0.8δ、1.0δ、1.5δ時Ex和Ey本文方法計算結果與解析解的振幅及誤差(圖8)。

根據CSAMT近區、過渡區和遠區的劃分原則可知，當加載范圍r=0.5δ時，場源完全處于近區范圍內。從圖8可以看出，當場源加載范圍r=0.5δ時，測站位于-320～320m范圍內，Ex和Ey數值解相對解析解出現明顯的跳變，最大相對誤差分別達到15%和10%，且其場源性質特征不明確，不能滿足計算精度的要求；當加載范圍r=0.8δ、測站位于-320～320m范圍內時，Ex和Ey數值解相對解析解也出現一定的跳變，但相對于r=0.5δ計算結果，Ex和Ey數值解的最大相對誤差分別降至5%和7%，且其場源性質表現為一定程度的連續變化的特征；當場源加載范圍r=1.0δ時，場源加載范圍已基本包含了近區范圍，Ex和Ey的數值解與解析解非常接近，不存在跳變，在測線中垂線附近，其相對誤差曲線規律性更強，誤差更小。除緊靠中垂線的幾個測點外，其他測點的數值模擬結果相對于r=0.5δ和r=0.8δ兩種情形的計算結果更優。因此，當加載范圍為r=1.5δ時，Ex和Ey數值解精度更高，特別是測線中垂線附近測點的相對誤差得到了很好的控制。因此，利用等效場源進行CSAMT三維數值模擬時，對于場源的加載范圍最好不小于1.5倍趨膚深度。

圖8 場源不同加載范圍時Ex(左)和Ey(右)數值解和解析解(a)及相對誤差(b)

6 結論

本文從CSAMT正演模擬三維邊值問題出發，采用等效場源法處理場源的奇異性，并利用無限元代替傳統截斷邊界條件，通過有限元-無限元耦合法，利用散度校正和直接法求解方程組的策略，實現了小區域的三維CSAMT快速高精度正演模擬。

與均勻半空間模型解析解的對比、分析可見，基于等效場源的CSAMT三維無限元數值模擬程序數值解與解析解平均相對誤差均小于1%，精度較高，驗證了基于等效場源法的三維CSAMT無限元方法的正確性；通過相同區域無限元法與傳統有限元法結果對比、無限元法與傳統大區域有限元法計算效率對比以及多頻點計算結果的對比，可以看出，本文提出的基于等效場源的三維無限元法具有離散區域小、求解速度快、計算精度高等優點。

針對等效場源法中精確加載場源的問題，數值模擬結果表明，基于等效場源的CSAMT三維無限元數值模擬結果，場源等效模擬最佳范圍應不小于1.5倍趨膚深度，其計算結果即能滿足精度要求。