李小康,王 偉,劉 磊

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

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

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

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

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

1 方法

1.1 方程離散化

角頻率為ω的自由頻率空間中,導電介質的麥克斯韋方程為:

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

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

體積與電導率的乘積為:

以棱邊 (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進行研究。

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

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

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

1.2 多重網格法

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

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

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

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

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

其中

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

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

2 模型正演計算

2.1 垂直接觸帶模型

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

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

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

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

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

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

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

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

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

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

2.2 均勻介質中電流源模型

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

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

表4 時點電流源模型計算結果

表5 時點電源模型計算結果

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

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

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

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

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

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

4 結論與思考

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

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

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

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

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

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

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

[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.