李小康,王 偉,劉 磊

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

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

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

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

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

1 方法

1.1 方程離散化

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

ωμ0σ～- ¤E ×μ-1r¤ ×E=-ωμ0Js (1)式中,矢量 E(ω,x)表示相應(yīng)角頻率和位置上的電場(chǎ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值可忽略不計(jì)。在后面算例中,我們針對(duì)地球物理模型取:r=1,而理想化模型中取:r=0。

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

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

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

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

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

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

經(jīng)過(guò)這樣的離散化,最終的解對(duì)于常系數(shù)具有二階精度[5],即便系數(shù)在跨域過(guò)程中發(fā)生變化,精度仍然能夠保持,除非發(fā)生的變化不連續(xù)。因?yàn)樵谏鲜銮闆r下,一般會(huì)導(dǎo)致解中存在一階的誤差。

1.2 多重網(wǎng)格法

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

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

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

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

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

其中

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

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

2 模型正演計(jì)算

2.1 垂直接觸帶模型

基于本征函數(shù)的垂直接觸帶問(wèn)題是經(jīng)常用于算法測(cè)試的模型,此處對(duì)其稍作修改,使用理想導(dǎo)電邊界條件。設(shè)研究域?yàn)?Ω= [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í)設(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ù),則電場(chǎng)的切向分量在邊界處為零。電流源表達(dá)式可寫(xiě)為 Ie=-δE+¤× (ωμ)-1¤×E,將精確解式 (3)代入該式得到:

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

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

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

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

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

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

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

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

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

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

表4 時(shí)點(diǎn)電流源模型計(jì)算結(jié)果

表5 時(shí)點(diǎn)電源模型計(jì)算結(jié)果

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

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

進(jìn)一步地為了構(gòu)建一個(gè)更加貼近實(shí)際的地下半空間模型,對(duì)已有的SEG/EA GE的鹽丘模型[8]進(jìn)行一定修改。修改前的鹽丘模型,是設(shè)計(jì)用來(lái)模擬人工地震波傳播,其中構(gòu)建了一個(gè)被沉積巖所包圍的復(fù)雜鹽丘體,設(shè)計(jì)海水深度120m,模型尺度設(shè)置為13500m×13480m×4680m。由于在實(shí)際的海洋可控源電磁勘探 (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ù)填入通過(guò)雙曲余弦拉伸的新網(wǎng)格。其中水平方向的網(wǎng)格剖分與原模型相同,設(shè)定z范圍為 (-5000m,5000m)。在海底,在y=0m處沿 x方向由 x=4000m至x=10000m之間,按照200m間距設(shè)置接收機(jī),深度范圍設(shè)定在605m至625m,模型的示意圖見(jiàn)圖4。為對(duì)模型進(jìn)行離散化,需計(jì)算相應(yīng)的元胞均值,頻率為

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

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

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

4 結(jié)論與思考

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

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

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

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

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

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

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

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