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

張林成 胡宏伶 湯井田 肖衛(wèi)初 肖 曉 原 源

( ①湖南城市學(xué)院信息與電子工程學(xué)院，湖南益陽 413002；②中南大學(xué)地球科學(xué)與信息物理學(xué)院湖南長沙 410083；③計(jì)算與隨機(jī)數(shù)學(xué)教育部重點(diǎn)實(shí)驗(yàn)室(湖南師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院)，湖南長沙 410081；④核資源與環(huán)境國家重點(diǎn)實(shí)驗(yàn)室(東華理工大學(xué))，江西南昌 330013 )

0 引言

隨著深部資源勘探的不斷深入，要求可探測深度越來越大[1]?？煽卦匆纛l大地電磁測深法 (CSAMT，Controlled Source Audio Electromagne-tic Method)[2]以勘探深度大、抗干擾能力強(qiáng)、工作效率高等特點(diǎn)成為研究熱點(diǎn)，是中深部資源勘探的重要手段之一。該方法在隱伏礦體勘察、復(fù)雜地區(qū)頁巖氣勘探、地?zé)豳Y源調(diào)查及海洋油氣勘探等領(lǐng)域都取得了良好效果[3-4]，成為發(fā)展快且前景可期的地球物理勘探方法之一。

CSAMT電磁場正演的主流方法主要有邊界單元法、有限差分法、積分方程法和有限單元法四種，其中有限單元法以其理論系統(tǒng)化、適應(yīng)性強(qiáng)、計(jì)算精度高、弱解可微等優(yōu)點(diǎn)得到了更多的重視和應(yīng)用。Coggon[5]于1971年提出了大地電磁問題的有限元模擬方法，自此有限元計(jì)算在電磁領(lǐng)域得到了極大的發(fā)展。CSAMT有限元數(shù)值模擬中，邊值問題包括控制方程和邊界條件兩個(gè)方面。對(duì)于控制方程，由于場源存在奇異性，場源的處理方式是關(guān)鍵，常用的方法有二次場法和總場法，其中二次場法是主流。二次場法將場分解為背景場和異常場，背景場利用均勻半空間或?qū)訝钅Ｐ徒馕鼋饪芍苯佑?jì)算，二次場則通過有限元法求取[6-12]?？倛龇ㄖ苯訌目倛鲋郑捎媒品M奇異性場源特征(例如偽delta函數(shù)法)，然后通過有限元求解場值[13-22]。在CSAMT三維有限元正演模擬中，無論采用總場法還是二次場法，場值的求解精度都是正演模擬是否成功的標(biāo)志，因此開展不同場源方案下場值求解精度的研究十分必要。

一般來說，對(duì)于邊界條件的處理基本上采用傳統(tǒng)截?cái)噙吔绶椒╗6-22]，即在一個(gè)相對(duì)較大的區(qū)域內(nèi)，將無窮邊界問題近似為有限區(qū)域，這往往會(huì)造成有限元計(jì)算區(qū)域太大、節(jié)點(diǎn)數(shù)太多、存儲(chǔ)量過大和計(jì)算耗時(shí)過長等問題。較合理的邊界處理策略應(yīng)該是將計(jì)算區(qū)域限制在較小的目標(biāo)區(qū)域內(nèi)，以減少單元數(shù)和節(jié)點(diǎn)數(shù)，實(shí)現(xiàn)節(jié)約計(jì)算資源及提高計(jì)算速度的目的。有限元-無限元耦合法是針對(duì)常規(guī)截?cái)噙吔缫鸬拇蟪叨饶Ｐ蛦栴}提出的[23]，具有離散區(qū)域小、求解速度快和計(jì)算精度高等優(yōu)點(diǎn)。其核心思路是用“無限單元”取代傳統(tǒng)截?cái)噙吔?，通過坐標(biāo)映射將“無限單元”沿某一方向無限擴(kuò)展，實(shí)現(xiàn)電磁場在無限單元中快速衰減為零，最后通過有限元-無限元耦合求解方程組，從而實(shí)現(xiàn)CSAMT的三維正演計(jì)算。

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

針對(duì)以上研究現(xiàn)狀，本文從基于總場方案的邊值問題出發(fā)，首先采用等效場源法，假設(shè)場源周圍很小區(qū)域內(nèi)電磁場的分布是均勻的，用區(qū)域內(nèi)網(wǎng)格節(jié)點(diǎn)上電磁場的近場值表示場源，并將節(jié)點(diǎn)場值作為第一類邊界條件帶入場方程，實(shí)現(xiàn)場源的近似模擬；然后，用“無限單元”代替?zhèn)鹘y(tǒng)截?cái)噙吔鐥l件，采用Garlerkin有限元法和Pardiso并行直接求解器實(shí)現(xiàn)了三維CSAMT快速正演計(jì)算；最后，通過均勻半空間模型解析解驗(yàn)證方法的正確性，并開展了場源等效模擬最佳范圍的研究。

1 邊值問題

1.1 基本方程

對(duì)于角頻率為ω、時(shí)間因子為e-iω的水平電偶極子，在各向同性導(dǎo)電介質(zhì)中產(chǎn)生的電場E滿足雙旋度方程[29]

(1)

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

1.2 場源的計(jì)算

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

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

在柱坐標(biāo)系中，上、下半空間中任意一點(diǎn)的電場分量為

(2)

(3)

根據(jù)電法理論，在電性分界面上，電場的法向分量是不連續(xù)的，而節(jié)點(diǎn)有限元要求電磁場法在法向上是連續(xù)的。因此，求解得到的有限元解往往不準(zhǔn)確，需要對(duì)解進(jìn)行校正。

在有源區(qū)，節(jié)點(diǎn)有限元的電場解不滿足

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

在無源區(qū)域不滿足

·(εE)=0

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

(4)

圖1 均勻半空間參數(shù)及坐標(biāo)示意圖

1.3 Galerkin有限元法

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

r=××E+k2E+·E

(5)

取任意測試函數(shù)V，在區(qū)域Ω上令

?Ωr·VdΩ=0

有

(6)

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

E=0

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

(7)

利用有限元對(duì)內(nèi)部計(jì)算區(qū)域進(jìn)行離散。假設(shè)內(nèi)部計(jì)算區(qū)域內(nèi)有n個(gè)節(jié)點(diǎn)，對(duì)于第j個(gè)測試函數(shù)Vj有

(8)

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

有限元—無限元耦合算法即將整個(gè)求解區(qū)域分為有限元區(qū)域和無限元區(qū)域，用無限元區(qū)域代替?zhèn)鹘y(tǒng)的外邊界，在這兩個(gè)區(qū)域分別采取有限元分析和無限元分析，通過總體剛度矩陣組裝將兩者結(jié)合，從而進(jìn)行數(shù)值求解。

圖2為有限元和無限元計(jì)算區(qū)域劃分示意圖。圖中灰黑色區(qū)域?yàn)橛邢拊獏^(qū)域，即目標(biāo)區(qū)域，包含場源、目標(biāo)體及測點(diǎn)等；白色區(qū)域?yàn)闊o限元區(qū)域，即從有限元區(qū)域邊界至無窮遠(yuǎn)，是邊界計(jì)算區(qū)域。無限元分析就是在該區(qū)域的某一方向上采用無限元映射和形函數(shù)，將整體坐標(biāo)映射到局部坐標(biāo)上，其原理與有限元分析相同。

圖2 有限元和無限元計(jì)算區(qū)域劃分示意圖

圖中P為坐標(biāo)原點(diǎn)，節(jié)點(diǎn)1～8為無限元基本要素，節(jié)點(diǎn)1～4為有限元邊界上某單元上的節(jié)點(diǎn)，P點(diǎn)到節(jié)點(diǎn)5～8的距離分別為P點(diǎn)到節(jié)點(diǎn)1～4距離的兩倍，圖4同

2.1 有限元分析

有限元單元分析時(shí)，采用矩形六面體進(jìn)行區(qū)域離散，單元節(jié)點(diǎn)編號(hào)及坐標(biāo)如圖3所示。

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

圖3中，子、母單元坐標(biāo)對(duì)應(yīng)關(guān)系為

(9)

矩形六面體形函數(shù)表達(dá)式為

(10)

式中(ξq，ηq，ζq)(q=1，2，…，8)表示子單元節(jié)點(diǎn)q在母單元中的坐標(biāo)。

2.2 無限元分析

無限元單元分析時(shí)，采用三維八節(jié)點(diǎn)Astley型無限元。圖4是三維無限元映射示意圖。無限元分析就是通過無限元映射將無窮坐標(biāo)映射到圖4b中的局部坐標(biāo)上。圖4b中無限單元最外圍的四個(gè)節(jié)點(diǎn)9～12代表無窮遠(yuǎn)，其場值為零。

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

圖4b中ξ方向表示無限元映射方向，在ζ-η平面內(nèi)，無限單元與有限單元采用相同的映射形式和形函數(shù)。無限元坐標(biāo)映射為

(11)

再結(jié)合ξ方向的坐標(biāo)映射關(guān)系式，可得到8個(gè)結(jié)點(diǎn)無限單元坐標(biāo)映射函數(shù)

(12)

無限元形函數(shù)M的表達(dá)式為

(13)

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

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

3 方程組求解

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

Ax=b

(14)

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

4 精度驗(yàn)證

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

4.1 模型計(jì)算

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

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

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

4.2 相同區(qū)域無限元法與傳統(tǒng)有限元法結(jié)果對(duì)比

針對(duì)上述均勻半空間模型，對(duì)相同區(qū)域進(jìn)行相同的網(wǎng)格剖分，利用傳統(tǒng)有限法進(jìn)行模型正演，并對(duì)比分析本文方法與傳統(tǒng)有限元法的計(jì)算結(jié)果。圖6為無限元法與傳統(tǒng)有限元法電場分量Ex和Ey振幅數(shù)值解與解析解的相對(duì)誤差。

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

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

4.3 無限元法與傳統(tǒng)大區(qū)域有限元法計(jì)算效率對(duì)比

關(guān)于無限元法與傳統(tǒng)大區(qū)域有限元法計(jì)算效率對(duì)比，文獻(xiàn)[10]有詳細(xì)分析，此處不再贅述。從表1可以看出，相對(duì)于基于總場法的傳統(tǒng)大區(qū)域有限元法，本文提出的基于等效場源的CSAMT三維無限元計(jì)算性能更高，明顯具有離散區(qū)域小、計(jì)算節(jié)點(diǎn)少、求解速度快等優(yōu)勢。

圖6 無限元法與傳統(tǒng)有限元法Ex(a)和Ey(b)數(shù)值解相對(duì)誤差

表1 基于總場法的傳統(tǒng)大區(qū)域有限元法與小區(qū)域無限元法計(jì)算數(shù)據(jù)對(duì)比[10]

4.4 無限元法多頻點(diǎn)計(jì)算結(jié)果對(duì)比分析

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

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

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

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

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

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

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

6 結(jié)論

本文從CSAMT正演模擬三維邊值問題出發(fā)，采用等效場源法處理場源的奇異性，并利用無限元代替?zhèn)鹘y(tǒng)截?cái)噙吔鐥l件，通過有限元-無限元耦合法，利用散度校正和直接法求解方程組的策略，實(shí)現(xiàn)了小區(qū)域的三維CSAMT快速高精度正演模擬。

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

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