多重網(wǎng)格算法在地球物理電磁三維正演中的應(yīng)用

李小康,王 偉,劉 磊

(1.中國地質(zhì)大學(xué) (北京)地球物理與信息技術(shù)學(xué)院,北京100083;2.中國地質(zhì)大學(xué) (北京)地球科學(xué)與資源學(xué)院,北京100083)

大地中電磁波的傳播可以用麥克斯韋方程組和歐姆定律進(jìn)行表述。在地球物理勘探一般用到的是1Hz左右的低頻部分,光速波的幅值很小,并且波長大大超出研究域范圍,因此擴(kuò)散作用是主要的研究對象,大地電磁法和可控源電磁法應(yīng)用的都是電磁波的擴(kuò)散原理。

在二維條件下,電磁波傳播問題可被簡化為泊松型方程,方程可以使用直接法求解,如嵌套剖分[1];也可以使用迭代法求解,如多重網(wǎng)格法。而在三維條件下,直接求解法對于大尺度問題進(jìn)行求解的成本太高,而迭代法就顯示出其優(yōu)勢。

在求解線性方程組時,旋度-旋度運算所形成的大型零空間是一個很大的難題。為了避免求解零空間問題,一般通過亥姆霍茲分解將電場轉(zhuǎn)為電位,得到泊松型方程組系。實際上,零空間的問題可以直接通過散度校正求解,或者通過劃分小的局部子系統(tǒng)進(jìn)行間接求解。本文選用后一種方法,用有限積分法對方程組使進(jìn)行離散化,這一過程可以被看作在 Yee方案的基礎(chǔ)上[2],引入了有限體積法。

多重網(wǎng)格法應(yīng)用于三維電磁法正演問題,一般有兩種思路,即單獨使用多重網(wǎng)格法求解,或者是作為 Krylov子空間法的預(yù)處理過程使用。當(dāng)多重網(wǎng)格法本身很難消除某種誤差的情況下,一般使用后一種方案。本文對于以上兩種方案的收斂性都進(jìn)行了研究,并進(jìn)行了一系列模型試算。進(jìn)行正演試算的模型包括,理想狀態(tài)下的垂直接觸帶模型,均勻介質(zhì)中點電源模型和海洋環(huán)境鹽丘模型。其中,對于網(wǎng)格拉伸對于算法收斂性的影響和解決方案,本文也進(jìn)行了著重研究。

1 方法

1.1 方程離散化

角頻率為ω的自由頻率空間中,導(dǎo)電介質(zhì)的麥克斯韋方程為:

ωμ0σ～- ¤E ×μ-1r¤ ×E=-ωμ0Js (1)式中,矢量 E(ω,x)表示相應(yīng)角頻率和位置上的電場分量,Js(ω,x)為電流源,σ～(X) =σωμrμ0,σ(X)為電導(dǎo)率,:r(X)為介電常數(shù),μr(X)為磁導(dǎo)率,:0和μ0為真空中的值,使用SI單位制。由于地球物理勘探中所應(yīng)用的頻率在1Hz左右,那么除在空氣中外,σ～(X)中的ω:0:r值可忽略不計。在后面算例中,我們針對地球物理模型取:r=1,而理想化模型中取:r=0。

下一步,在張量積笛卡爾網(wǎng)格中,使用有限積分法對方程組進(jìn)行離散化。網(wǎng)格節(jié)點可以確定方形網(wǎng)格,電場分量取棱邊上的均值,而磁場分量 H= (¤ ×E)/(ωμrμ0)通常被分配指向元胞的各個表面。這些分量可以形成節(jié)點位于原始元胞網(wǎng)格質(zhì)心的派生網(wǎng)格,見圖1。磁場分量的旋度可以被歸到原始元胞的各個棱邊,但原始元胞的各節(jié)點并不一定要位于派生網(wǎng)格的質(zhì)心。

體積與電導(dǎo)率的乘積為:

以棱邊 (k+1/2,l+1/2,m+1/2)處的分量為例,求平均過程為

Vk+1/2,l,m是兩個派生元胞疊加的部分,其中包含了棱邊 (k+1/2,l,m),也是四個原始元胞的公共邊。為了運算方便,我們針對 Sk+1/2,l,mE1k+1/2,l,m進(jìn)行研究。

由電場旋度求取電場分量,需要求元胞各表面處的 1/μr,即在派生元胞棱邊處對 1/μr求平均[4]。由于物性參數(shù)值在每一個元胞中為常數(shù),這個過程也就是針對平面兩側(cè)元胞中的1/μr按體積求加權(quán)平均。有時為了使模擬更加精確的逼近于真實情況,往往設(shè)定元胞內(nèi)部的物性也存在變化,取平均的過程就不能使用簡單的體積加權(quán)。這樣的情況較復(fù)雜,本文暫時不進(jìn)行考慮。

對于電流源Js的各個分量可以按照與電場分量類似的方法分配到各個棱邊上,各個分量大小分別為派生元胞各個表面的均值,磁場分量的旋度處理方式與之相同。如果將Js通過派生元胞一個表面的分量的積分與其所在的棱邊長度相乘,那么就得到了電流源分量與雙重棱邊體積的離散化結(jié)果。

經(jīng)過這樣的離散化,最終的解對于常系數(shù)具有二階精度[5],即便系數(shù)在跨域過程中發(fā)生變化,精度仍然能夠保持,除非發(fā)生的變化不連續(xù)。因為在上述情況下,一般會導(dǎo)致解中存在一階的誤差。

1.2 多重網(wǎng)格法

Feigh曾介紹各類多重網(wǎng)格方法[6],本文網(wǎng)格粗化使用八個相鄰細(xì)網(wǎng)格元胞自然粗化為一個粗網(wǎng)格元胞的方法,引入將物性參數(shù)與相應(yīng)元胞體積的乘積用于約束。為了約束殘差,共引入兩個約束算子。

第一個等于在粗網(wǎng)格雙重棱邊體積內(nèi)包含的細(xì)網(wǎng)格雙重棱邊體積部分的縮放比例合適情況下,給定棱邊處的粗網(wǎng)格體積中包含的所有雙重棱邊體積的積分。如果元胞的寬度固定,那么這個值

圖1 原始元胞 (空心)與派生元胞 (實心)示意圖

物性參數(shù)在每個元胞中為常數(shù),要得到電場分量所在各邊上的電導(dǎo)率σ,需取均值。在兩個派生元胞共享的平面上,應(yīng)當(dāng)使用面積平均[3],下面詳述取均值的過程。

將一個 Nx×Ny×Nz的網(wǎng)格的節(jié)點分別指定為xk(k=0,…,Nx),yl(l=0,…,Ny)和zm(m=0,…,Nz)。若某個元胞的體積表示為:則減少到全權(quán)重。為了對其進(jìn)行詳細(xì)闡述,首先考慮一種較為簡單的情況,有一個 Nx×Ny×Nz的網(wǎng)格,其中 Nx=Ny=Nz=2M,M≥1。假定按照網(wǎng)格節(jié)點 xk(k=0,…,Nx),yl(l=0,…,Ny)和zm(m=0,…,Nz)中隔點對網(wǎng)格進(jìn)行粗化,粗化后的網(wǎng)格節(jié)點為 (x2K,y2L,z2M)其中K=0,…,Nx/2,L=0,…,Ny/2,M=0,…,Nz/2。以殘差 r的第一分量為例,表述約束因子。將殘差與雙重棱邊體積相乘,其分量與相應(yīng)的電場分量分布在相同的棱邊上。第一級粗網(wǎng)格的殘差分量就是,上標(biāo)2h指示網(wǎng)格的粗化值。令粗網(wǎng)格中的指標(biāo) K與細(xì)網(wǎng)格中存在k=2 K的關(guān)系,同樣的l=2L,m=2M。那么該分量的約束條件為:

其中

第二個約束因子,等于在雙重棱邊體積縮放比例合適的情況下,粗網(wǎng)格棱邊上細(xì)網(wǎng)格殘差的

為了將粗化網(wǎng)格時電場分量的校正值進(jìn)行延拓,在分量方向上分段使用常數(shù)插值,在垂直于分量的平面上使用線性與雙線性插值,這也與式(2)的共軛方程相同。在編序過程中使用對稱體Gauss-Seidel(SBGS)法,作為網(wǎng)格循環(huán)過程的平滑器和SBGS松弛的一個步驟,每個循環(huán)進(jìn)行兩次后平滑處理。而當(dāng)多重網(wǎng)格充當(dāng)預(yù)處理器時,僅需作一次網(wǎng)格循環(huán)。

2 模型正演計算

2.1 垂直接觸帶模型

基于本征函數(shù)的垂直接觸帶問題是經(jīng)常用于算法測試的模型,此處對其稍作修改,使用理想導(dǎo)電邊界條件。設(shè)研究域為,Ω= [0.2π]3m3,ψ=sinkxsinlysimy,其中 k、l、m為正整數(shù)。精確解為:

E1=α19xψ,E2= α29yψ,E3=α39zψ (3)

設(shè)研究域Ω中有,Ω1(z<π)和Ω1(z>π)。在Ω1中σ=σ0+σ1(x+1) (y+2) (z-π)2,在Ω2中σ=σ0,同時設(shè):r=0,μr=1,ω=106H z。其他參數(shù)設(shè)置為α1=α2=-2V,α2=1V,k=l=m=1,σ0=105/m,σ1=15/m。由于模型中使用理想導(dǎo)體邊界條件,基于正弦函數(shù),則電場的切向分量在邊界處為零。電流源表達(dá)式可寫為 Ie=-δE+¤× (ωμ)-1¤×E,將精確解式 (3)代入該式得到:

同時,將三個方向上網(wǎng)格剖分單元數(shù)量分別設(shè)為 Nx=Ny=Nz=N,對于每套網(wǎng)格分別使用兩種迭代方法,即單獨使用多重網(wǎng)格迭代 (M G)和多重網(wǎng)格預(yù)處理后的雙共軛梯度穩(wěn)定法 (bi),設(shè)定當(dāng)殘差對于零解的范數(shù)由初始值下降到10-8數(shù)量級時,停止迭代,結(jié)果見表1。

通過2范數(shù)和最大范數(shù)對誤差進(jìn)行檢查。如果網(wǎng)格有 Nx×Ny×Nz個元胞,網(wǎng)格節(jié)點可表示為(xk,yl,zm),k=0, …,Nx,l=0, …,Ny,m=0,…,Nz,第一電場分量的2范數(shù)為:其他分量的范數(shù)與之類似。解的誤差為:相對誤差為:最大范數(shù)為:

表1 垂直接觸帶模型迭代次數(shù)及誤差統(tǒng)計

由表1可見,此時多重網(wǎng)格法的收斂性與網(wǎng)格尺度無關(guān)。雙共軛梯度穩(wěn)定法的迭代次數(shù)較少,但是由于額外的CPU時間和存儲要求,使得該方法變得十分不合算。例如,對于1283尺寸的網(wǎng)格而言,雙共軛梯度穩(wěn)定法6次迭代與單獨使用多重網(wǎng)格算法8次迭代相比,效率僅高10%。通過對誤差的分析,確定了數(shù)值解的二階精度。

接下來對網(wǎng)格拉伸的效果進(jìn)行研究,網(wǎng)格拉伸是指從坐標(biāo)原點沿坐標(biāo)軸開始向外進(jìn)行網(wǎng)格剖分時,按照1+α的比率將剖分尺度逐步放大的處理,等距網(wǎng)格中α=0。針對拉伸系數(shù)為α=0.04的網(wǎng)格剖分模型計算結(jié)果見表2。通過對比,顯然可以發(fā)現(xiàn)此時多重網(wǎng)格法收斂性對于網(wǎng)格剖分的獨立性消失,而雙共軛梯度穩(wěn)定法的優(yōu)勢也顯現(xiàn)出來,因為該方法可以對收斂較慢的分量進(jìn)行有效地處理。由表中誤差檢測的結(jié)果:2和:max可知,兩種方法的解仍具備二階精度。

表2 α=0.04時垂直接觸帶模型計算結(jié)果

綜上可知,對于電流源問題,網(wǎng)格拉伸未必會起到好的效果,至少一定程度上會影響迭代過程的收斂速度,在未來的研究過程中應(yīng)當(dāng)引起注意。

2.2 均勻介質(zhì)中電流源模型

下面對更加貼近實際的勻質(zhì)地層中的點電流源問題進(jìn)行討論。模型參數(shù)設(shè)置為研究域設(shè)定成尺寸為[-1,1]×[-1,1]×[-1,1]km3的立方體區(qū)域,點J3= (0.01)δ (x)A/m2電流源放置在原點。首先采取等距網(wǎng)格對模型進(jìn)行離散化,剖分單元設(shè)置為Nx=Ny=Nz=N,然后使用功率拉伸網(wǎng)格對兩種迭代法的表現(xiàn)進(jìn)行對比。為了電流源表達(dá)式進(jìn)行離散化,這里使用了三線性插值共軛法。迭代停止條件是數(shù)值解對于精確解殘差的l2范數(shù)由初始值降至10-8數(shù)量級,該問題的精確解由Ward給出[7]。盡管理想導(dǎo)體邊界條件不適用于該問題,但是只要將邊界置于據(jù)電流源足夠遠(yuǎn)處,在求解時仍然可以使用。由于精確解在源處存在一個r-3的奇異值,所以無法通過網(wǎng)格加密降低誤差。點電流源可以由有限長線源代替,線源的長度則取決于源所在的棱邊長度。在本算例中,電源的位置恰好位于節(jié)點處,并恰好可以由節(jié)點上方和下方的棱邊表示。上述棱邊長度的改變對近區(qū)的電磁場影響較大,有可能使最終的解失去二階精度。而遠(yuǎn)區(qū)場則不受源的這些細(xì)節(jié)影響,只取決于源的總體性質(zhì)。通過測量立方體 [-250,250]3m3之外的誤差,得到的計算結(jié)果見表3～表5。

表3 點電流源模型固定網(wǎng)格的迭代次數(shù)及誤差統(tǒng)計

表4 時點電流源模型計算結(jié)果

表5 時點電源模型計算結(jié)果

可見在網(wǎng)格無拉伸的情況下,迭代效率相當(dāng)高,而對網(wǎng)格進(jìn)行超過幾個百分點的拉伸后,收斂性就開始急劇下降。

2.3 海洋環(huán)境鹽丘模型

進(jìn)一步地為了構(gòu)建一個更加貼近實際的地下半空間模型,對已有的SEG/EA GE的鹽丘模型[8]進(jìn)行一定修改。修改前的鹽丘模型,是設(shè)計用來模擬人工地震波傳播,其中構(gòu)建了一個被沉積巖所包圍的復(fù)雜鹽丘體,設(shè)計海水深度120m,模型尺度設(shè)置為13500m×13480m×4680m。由于在實際的海洋可控源電磁勘探 (CSEM)中,海水深度比原模型設(shè)置要深,因此首先需要把模型中海水深度加大到500m。然后由于需要進(jìn)行電磁三維正演,因此物性參數(shù)中的速度相應(yīng)由電導(dǎo)率代替,即將海水中波速1500m/s替換為電導(dǎo)率10/3 S/m,鹽丘中波速4000m/s替換為電導(dǎo)率1/30 S/m,而深度大于3956m處的基底巖層電導(dǎo)率設(shè)為500 S/m,沉積巖中的電導(dǎo)率由波速 v(單位為m/s)確定[9],設(shè)為空氣中的電導(dǎo)率設(shè)為10-10S/m,最后將這些電導(dǎo)率參數(shù)填入通過雙曲余弦拉伸的新網(wǎng)格。其中水平方向的網(wǎng)格剖分與原模型相同,設(shè)定z范圍為 (-5000m,5000m)。在海底,在y=0m處沿 x方向由 x=4000m至x=10000m之間,按照200m間距設(shè)置接收機(jī),深度范圍設(shè)定在605m至625m,模型的示意圖見圖4。為對模型進(jìn)行離散化,需計算相應(yīng)的元胞均值,頻率為

源位于 (x,y,z)= (6400m,6500m,500m)與 (x,y,z)= (6600m,6500m,500m)之間。

針對1283剖分網(wǎng)格進(jìn)行計算的結(jié)果見表6,網(wǎng)格拉伸由雙曲余弦函數(shù)控制,其中最小的單元格寬度為hmin,相鄰兩網(wǎng)格之間的尺度比的最大值用表示。正演計算過程中,多重網(wǎng)格還是作為雙共軛梯度穩(wěn)定法的預(yù)處理步驟使用,同樣地,網(wǎng)格拉伸對于迭代性能的影響也做了進(jìn)一步研究。

表6 使用多重網(wǎng)格預(yù)處理的bicgstab法迭代次數(shù)

4 結(jié)論與思考

本文研究了針對在張量積笛卡爾網(wǎng)格中的電磁擴(kuò)散方程有限積分離散化技術(shù)的多重網(wǎng)格法,進(jìn)行了一系列模型試算,著重觀察了等距網(wǎng)格和拉伸網(wǎng)格剖分對收斂性的影響。網(wǎng)格拉伸會引起離散方程的各向異性的特征,這一直以來都是多重網(wǎng)格法比較難以解決的問題,未來的研究中可以考慮以下幾種解決方案:

1)Krylov子空間法可以用于消除導(dǎo)致收斂變慢的解分量。盡管本文使用雙共軛梯度穩(wěn)定法取得了一定效果,但是該方法起到的作用有限。

2)線性松弛法,例如對稱線性 Gauss-Seidel法,可以用來消除網(wǎng)格拉伸帶引起的各向異性趨向。對于三維問題,可以考慮使用平面松弛法,但是計算成本會提高。

3)半粗化處理也可以達(dá)到同樣的效果。在標(biāo)準(zhǔn)的多重網(wǎng)格迭代流程中,每次僅對一個方向進(jìn)行粗化。使用對網(wǎng)格所有方向同時進(jìn)行粗化的算法,那么與之相關(guān)的時間成本增加可以通過后續(xù)計算效率的提高得以彌補(bǔ)。

4)通過權(quán)重約束因子和延拓因子的控制,可以大大降低物性參數(shù)變化帶來的負(fù)面效果[10]。

5)可以借鑒計算機(jī)技術(shù)中解決等距網(wǎng)格中局部細(xì)化問題的方法。這一算法應(yīng)當(dāng)可以在局部使用較細(xì)的網(wǎng)格,同時保證算法的收斂速度,當(dāng)然這一算法的編程實現(xiàn)更為復(fù)雜。

以上提到的幾種思路,必然會使得多重網(wǎng)格迭代的計算過程更為復(fù)雜,下一步需要做的就是從中選出效率或者說“性價比”最高的一種解決方案。

[1] George A,Liu J.Computer Solution of Large Sparse Positive Definite Systems[M]. Prentice-Hall,Englewood Cliffs,1981.

[2] Yee K.Numerical solution of initial boundary value p roblems involving Maxwell’s equations in isotropic media[J].Antennas Propag,1966,16:302-307.

[3] Clemens M.,Weiland T.Discrete electromagnetism with the Finite Integration Technique[J].Prog Electromagn Res(PIER),2001,32:65-87.

[4] van Rienen U.Frequency domain analysisof waveguidesand resonators with FIT on non-orthogonal triangular grids[J].Prog Electromagn Res(PIER) ,2001,32:357-381.

[5] Monk P,Süli E.A convergence analysisof Yee’s schemeon non uniform grids[J].SIAM J,1994,31:393-412.

[6] Butler DK.Potential Fields Methods for Location of Unexploded Ordnance[J].The Leading Edge,2001,20(8):890-8951.

[7] Ward S,H Hohmann.Electromagnetic theory for geophysical applications.Society of Exploration Geophysicists,Tulsa.Nabighian,M.N.(ed.)Electromagnetic Methods in Applied Geophysics—Theory,1987,1(1):131-311.

[8] Aminzadeh F,Brac J,Kunz T.3-D Salt and over thrust models[J]. Society of Exploration Geophysicists,Tulsa,1997.

[9] Meju M.A.,Gallardo L.A.,Mohamed A.K.Evidence for correlation of electrical resistivity and seismic velocity in heterogeneous near-surface materials[J].Geophys.Res.Lett.,2003,30(7):1373-1376.

[10] Wesseling P.An Introduction to Multi-grid Methods[M].Wiley,New York,1992.