周建美 劉文韜 魯凱亮 李 貅

(長安大學(xué)地質(zhì)工程與測繪學(xué)院，陜西西安 710054)

0 引言

大地電磁(MT)方法廣泛應(yīng)用于地殼研究[1-2]、巖石圈地幔調(diào)查[3-6]、資源勘探[7-11]和海洋調(diào)查[12]等領(lǐng)域。MT方法是一種典型的頻域電磁測深方法[13]，工作頻率范圍可達(dá)10-5～104Hz[14]。MT三維正演方法可以分為以下幾類：積分方程法[15-17]，有限差分方法[18-23]，有限體積方法[24-25]，有限元方法[26-32]以及無網(wǎng)格法[33]等。無論采用哪種正演方法，都需要求解多個頻率的大型稀疏線性方程組[34]。MT響應(yīng)頻率范圍通常很大，可能達(dá)到5個或更高數(shù)量級，每個數(shù)量級范圍內(nèi)通常包含5～10個采樣頻點(diǎn)[35]。此外，在每個頻率上都需要針對x方向TE極化源和y方向TM極化源進(jìn)行兩次正演[27]。因此，一次MT正演一般需要求解幾十個大型稀疏線性方程組，計(jì)算量很大。

學(xué)者們對此提出了多頻電磁正演加速方法。Um等[36]引入了一種針對迭代法的預(yù)處理重復(fù)利用技術(shù)，只需求解最低頻率的ILU預(yù)處理矩陣，并將該預(yù)處理矩陣重復(fù)應(yīng)用于其他頻率的正演計(jì)算，從而減少正演計(jì)算時間。但該算法對于較大頻率間隔情況下(比如MT)的正演模擬的有效性還有待驗(yàn)證。另一種常用方法是并行計(jì)算技術(shù)[37]，即采用多CPU計(jì)算設(shè)備并行計(jì)算多個頻率的正演響應(yīng)[38-39]。上述兩種改進(jìn)算法本質(zhì)上都是基于不同發(fā)射頻率正演響應(yīng)的獨(dú)立求解，當(dāng)發(fā)射頻點(diǎn)數(shù)增加時，方程求解次數(shù)或CPU數(shù)量也需要相應(yīng)增加。

相對上述兩種方法，Krylov子空間算法[40-43]更有效。隨頻率變化的電場響應(yīng)可以表示為一個傳遞函數(shù)與源矢量的乘積[40]，Krylov子空間算法可對其有效求解。Druskin等[40]基于多項(xiàng)式Krylov子空間近似的譜Lanczos方法(SLDM)求解上述電場響應(yīng)表達(dá)式，但該算法的收斂速度依賴于頻率范圍和模型的電導(dǎo)率差異[41]。相對而言，有理Krylov子空間方法[41-44]的收斂性更好，且收斂性不依賴于模型電導(dǎo)率差異和空間網(wǎng)格分布。但是該算法需要求解多個極點(diǎn)的大型線性方程組。采用重復(fù)極點(diǎn)，結(jié)合直接求解器，可以顯著提高求解速度。

盡管Krylov子空間算法可以求解人工源頻率域多頻電磁響應(yīng)，但是該方法不能直接推廣到MT正演求解。MT正演假設(shè)源是無限大的平面波源，通常采用邊界條件隱式代替源的作用[18-19]，相應(yīng)的MT正演響應(yīng)不能表示為一個傳遞函數(shù)與矢量的乘積形式。Kordy等[45]基于二次場方法，將MT響應(yīng)表示為一個傳遞函數(shù)與二次場矢量的乘積形式，從而可以采用Krylov子空間算法進(jìn)行求解，但是由于二次場是隨頻率變化的，因此求解每一個頻率的正演響應(yīng)時都需要重新構(gòu)建一次Krylov子空間。為此，Kordy等[45]引入一種有理插值技術(shù)解決這一問題，但仍需要多次正演計(jì)算生成有理插值方法所需的訓(xùn)練集，因而相比常規(guī)的逐次求解各個頻率正演響應(yīng)的方法，該方法僅僅獲得2～4倍的加速。

不同于Kordy等的處理方法，本文將MT的源顯式表示為一個位于高空的平面電流的形式[46]，該源對于所有頻率都是一個常數(shù)。當(dāng)平面電流位于模型上邊界時，與Zhang等[47]直接在模型上邊界賦值電場為1等效。在此基礎(chǔ)上，可以將MT響應(yīng)表示為一個傳遞函數(shù)與常數(shù)源矢量的乘積形式，從而可以采用有理Krylov方法求解MT響應(yīng)[47]。

與可控源電磁法正演問題不同[42]，模擬MT響應(yīng)需對每一個頻率計(jì)算TE和TM極化源的兩次正演，這等價(jià)于一個多源正演問題。塊狀Krylov技術(shù)廣泛應(yīng)用于多源正演問題，可有效提高計(jì)算效率[48-50]。為進(jìn)一步提高M(jìn)T正演的計(jì)算速度，不同于Zhang等[47]直接采用有理Krylov方法進(jìn)行求解，本文引入塊狀有理Krylov技術(shù)[51]，將TE和TM極化源表示為塊狀源矢量，將求解兩個極化源的正演響應(yīng)簡化為構(gòu)建一個塊狀有理Krylov子空間。

極點(diǎn)是有理Krylov方法的一個重要參數(shù)，極點(diǎn)的選取會影響方法的計(jì)算速度和計(jì)算精度。Güttel[52]對極點(diǎn)的選取算法進(jìn)行了總結(jié)，提出了一種最優(yōu)化極點(diǎn)選取方法，即通過全局搜索得到多個最優(yōu)化重復(fù)極點(diǎn)[50]?；谠摲椒?，本文引入一種更簡單、有效的方法獲取最優(yōu)化單個重復(fù)極點(diǎn)[42]，結(jié)合直接法求解器[53]，可實(shí)現(xiàn)MT三維模型的高效正演模擬。

本文首先介紹基于擬態(tài)有限體積的空間離散化方法[54]，源項(xiàng)顯式表示為面電流分布，從而將MT響應(yīng)表示為傳遞函數(shù)與常數(shù)源矢量乘積的形式； 然后，使用塊狀有理Krylov(BRK)算法獲得傳遞函數(shù)的模型降階近似解； 最后，通過典型的一維半空間模型和3D都柏林測試模型1(DMT1)正演[34]，驗(yàn)證本文算法的計(jì)算精度和計(jì)算速度。

1 正演算法

1.1 控制方程離散

設(shè)電磁場隨時間的變化關(guān)系為eiωt，忽略位移電流，則MT正演對應(yīng)的控制方程為[14]

?×E+iωB=0

(1)

?×μ-1B-σE=S

(2)

式中：E和B分別表示電場強(qiáng)度矢量和磁場強(qiáng)度矢量； 角頻率ω=2πf，其中f是頻率；μ是磁導(dǎo)率；σ是電導(dǎo)率；S是外加源項(xiàng)，顯式表示為平面電流分布[46]。采用自然邊界條件，利用基于交錯網(wǎng)格的擬態(tài)有限體積方法[54]，得到空間離散后的控制方程

Ce+iωb=0

(3)

CTMμb-Mσe=s

(4)

式中：C是離散形式的旋度算子；Mμ和Mσ分別是包含模型磁導(dǎo)率和電導(dǎo)率信息的離散矩陣，具體形式參見文獻(xiàn)[42]；e、b和s分別為離散形式的電場、磁場和發(fā)射源項(xiàng)。對于三維MT正演，同時考慮TE和TM極化源，采用如圖1所示的平面電流近似，平面電流位于空氣層中一定高度，可得空間離散的源項(xiàng)為

圖1 TE極化模式(左圖中藍(lán)色箭頭)和TM極化模式(右圖中橙色箭頭)的平面電流源分布示意圖

[H(x-xi)-H(x-xi+1)]

i=1,2,...,nx-1;j=1,2,...,ny

(5)

[H(y-yj)-H(y-yj+1)]

i=1,2,...,nx;j=1,2,...,ny-1

(6)

s=[s1,s2]

(7)

消去式(4)中的磁場，得到關(guān)于電場的控制方程

(A+iωI)e=-iωX

(8)

e(ω)=gω(A)X

(9)

1.2 塊狀有理Krylov方法

采用Krylov子空間方法近似求解gω(A)X，由于X是一個塊狀矩陣，因此可采用BRK方法一次性求解所有源的正演響應(yīng)[51]。給定塊狀矩陣X和單個重復(fù)極點(diǎn)ξ，BRK方法采用Gram-Schmidt正交化方法求解得到有理Krylov子空間的基，由一系列的塊狀矢量v1,v2,…,vm+1構(gòu)成，這里m表示子空間的維度。針對單個矢量和塊狀矩陣的有理Krylov方法(算法1)的程序?qū)崿F(xiàn)見表1[50]。

表1 算法1——塊狀有理Krylov方法

(10)

為了實(shí)現(xiàn)算法1，還需要預(yù)先給出極點(diǎn)ξ和迭代次數(shù)m。本文采用單個重復(fù)極點(diǎn)ξ，其值可通過近似收斂率公式快速得到[42]。對于頻率范圍T=[ωmin,ωmax]，最優(yōu)化極點(diǎn)為

(11)

得到最優(yōu)化極點(diǎn)ξ后，可以采用替代算法[42]快速得到最優(yōu)化迭代次數(shù)mopt。根據(jù)周建美等[42]的分析，頻率為ωmin和ωmax時，收斂率最差(模型降階解滿足給定殘差閾值所需迭代次數(shù)最大)，但是最低頻率ωmin對應(yīng)的相對誤差最大，因此可以簡單地選擇最低頻率ωmin對應(yīng)的模型降階解的相對誤差不再降低時的m作為mopt。由于已知最優(yōu)化的單個重復(fù)極點(diǎn)ξ，本文求解mopt的替代算法較B?rner等[41]提出的方法更加簡化，算法的具體計(jì)算機(jī)實(shí)現(xiàn)過程算法2如表2所示，其中：L0=diag(λ1,λ2…,λN)為替代對角矩陣，λ1,λ2，…,λN是[λmin=10-8,λmax=108]范圍內(nèi)N=500個對數(shù)等間隔分布的值；mmax為一個足夠大的數(shù)，本文取mmax=200；R(m)表示Rm的集合; arg minR(m)表示尋找R(m)中最小Rm所對應(yīng)的m。

表2 算法2——求解mopt的替代算法

圖2 fmax/fmin固定時不同fmin條件下有理Krylov解的相對誤差error隨m的變化曲線

圖3 求解頻率范圍的fmax固定時不同fmin的有理Krylov解的相對誤差error隨m的變化曲線

1.3 MT響應(yīng)

張量MT勘探典型的測量參數(shù)是地表水平電場Ex、Ey及全磁場分量Hx、Hy、Hz。假定有兩個正交極化源，阻抗張量為

(12)

其中下標(biāo)“1”和“2”分別代表TE和TM極化模式?；谧杩箯埩縕定義視電阻率和相位[13]分別為

(13)

(14)

式中：下標(biāo)“q”和“p”代表x或y； Re(·)和Im(·)分別表示取實(shí)部和虛部。

2 數(shù)值實(shí)例

通過兩個典型模型的正演驗(yàn)證BRK方法對MT正演計(jì)算的速度和精度。數(shù)值實(shí)例中空氣層電阻率均為1×106Ω·m。本文計(jì)算設(shè)備為32G內(nèi)存、四核主頻3.6GHz的Intel i7CPU的臺式電腦。

2.1 均勻半空間模型

采用BRK方法計(jì)算MT響應(yīng)。根據(jù)式(11)得到最優(yōu)化極點(diǎn)ξ=-0.1987。采用替代算法求解mopt。圖4為頻率fmin為0.001Hz時的有理Krylov子空間解的相對誤差error隨m的變化曲線。由圖可見，m>70時，error小于1×10-6且基本恒定，滿足計(jì)算精度要求，因此選取mopt=70。

圖4 均勻半空間模型fmin=0.001Hz的有理Krylov解的error曲線

圖5為不同h時計(jì)算的視電阻率和相位正演結(jié)果與解析解的對比。對于均勻半空間模型，TE和TM模式的響應(yīng)是相同的，因此圖5僅給出了TE模式的正演結(jié)果。

由圖5可見，在低頻時，靠近地表(h=0.5δ)時，平面電流源不能較好地近似為平面波源，因此正演響應(yīng)存在較大誤差； 當(dāng)平面電流源逐漸遠(yuǎn)離地表(h=4δ，8δ，12δ)時，正演誤差逐漸減??； 當(dāng)平面電流源位于空氣層邊界(h=16δ)時，所有頻點(diǎn)視電阻率相對誤差均小于1%，相位誤差均小于0.5°。

圖5 半空間模型不同h的MT視電阻率(a)和相位(b)響應(yīng)與解析解對比

均勻半空間模型計(jì)算結(jié)果表明，將平面電流源置于空氣層邊界時，本文BRK方法能夠精確地模擬MT三維正演響應(yīng)。

為了分析BRK方法的計(jì)算效率，對比本文BRK方法、非塊狀有理Krylov方法[42](RK)和逐次采用直接求解器PARDISO[53]求解各個頻率的正演響應(yīng)的直接求解法(Direct)的計(jì)算結(jié)果。雖然一維模型的TE模式和TM模式的響應(yīng)是相同的，但是為了通用性考慮，在比較算法計(jì)算效率時，本文依然以計(jì)算TE和TM模式響應(yīng)的總計(jì)算時間為依據(jù)。上述三種算法計(jì)算的MT響應(yīng)基本相同，因此本文僅對計(jì)算時間進(jìn)行比較，結(jié)果如表3所示?？梢?，本文BRK方法只需1次系數(shù)矩陣分解及70次矩陣回代； RK方法分別求解TE和TM模式的響應(yīng)，因此需要1次系數(shù)矩陣分解和140次矩陣回代； 直接求解法(Direct)需對每個頻率進(jìn)行1次系數(shù)矩陣分解，因此共計(jì)需要31次系數(shù)矩陣分解和31次矩陣回代。因此，本文BRK方法相比直接求解法，能夠顯著提高計(jì)算速度，加速比高達(dá)22.5； 相比RK方法，由于減小了1次子空間構(gòu)建的運(yùn)算次數(shù)，也能實(shí)現(xiàn)一定的加速效果。

表3 均勻半空間模型MT響應(yīng)模擬不同算法計(jì)算量及耗時對比

2.2 三維DTM1模型

為了進(jìn)一步分析BRK方法的精度和速度，利用三維都柏林測試模型1(DTM1)[30]進(jìn)行正演響應(yīng)計(jì)算，并與有限元法和積分方程方法的正演結(jié)果進(jìn)行對比。

DTM1模型的均勻半空間包含3個塊狀電性異常體，具體參數(shù)見表4和圖6。坐標(biāo)原點(diǎn)位于塊體1中心點(diǎn)的正上方地表，觀測點(diǎn)位于坐標(biāo)原點(diǎn)，fmin=10-4Hz，fmax=10Hz，在計(jì)算頻率范圍內(nèi)對數(shù)等間隔取21個頻點(diǎn)。計(jì)算區(qū)域離散為72×74×72個網(wǎng)格，最小網(wǎng)格長度為20m，網(wǎng)格放大系數(shù)為1.5。

圖6 DTM1模型示意圖

表4 DTM1模型塊體空間坐標(biāo)和電阻率

根據(jù)式(11)得到最優(yōu)化極點(diǎn)ξ=-0.1987，與一維模型的最優(yōu)化極點(diǎn)值相同，因?yàn)闃O點(diǎn)取決于最低頻率和最高頻率的乘積。采用替代算法求解mopt。圖7為fmin=1×10-4Hz的有理Krylov子空間解的相對誤差error隨m的變化曲線。相比一維模型，三維DTM1模型的求解頻率范圍更大，因此最小頻率fmin對應(yīng)的誤差收斂率更小，從而相對誤差穩(wěn)定時所對應(yīng)的m更大。由圖7可知，m>150時，相對誤差基本恒定，因此選取mopt=150。

圖7 fmin=1×10-4Hz時有理Krylov子空間解的error隨m的變化曲線

采用本文BRK方法得到的DTM1模型正演結(jié)果見圖8。文獻(xiàn)[32]給出了多種方法的DTM1模型正演響應(yīng)，其中積分方程法(IE)[55]和有限元法(FE)[27]的結(jié)果基本重合，因此本文以這兩種算法的正演結(jié)果作為參考，驗(yàn)證本文算法的精度。由圖8可知，BRK方法正演得到的xy和yx分量的視電阻率(ρxy，ρyx)和相位曲線(φxy，φyx)，與積分方程法[55]和有限元法[27]的結(jié)果基本重合； BRK方法正演得到的xx和yy分量的視電阻率(ρxx，ρyy)和相位(φxx，φyy)，在高頻段(>1Hz)由于視電阻率值太小沒有顯示，相位分布也相對雜亂； 但是在低頻段(<1Hz)，本文計(jì)算結(jié)果與積分方程法和有限元法的結(jié)果基本一致。因此，本文BRK方法能夠較精確地正演MT三維復(fù)雜模型的電磁響應(yīng)。

統(tǒng)計(jì)本文BRK法、RK[42]法和Direct法[53]開展DTM1模型正演的耗時如表5所示。相比Direct法，BRK法加速比達(dá)到11.8； 相比RK方法，也能實(shí)現(xiàn)一定的加速效果。相比一維半空間模型，三維DTM1模型的正演加速比有所降低，主要原因在于三維DTM1模型求解的頻率數(shù)量僅為21，每個數(shù)量級范圍內(nèi)僅選取4個頻率，相比一維半空間模型(頻率數(shù)量為31)有所減少，因而直接求解法所需的計(jì)算時間相對減少，導(dǎo)致BRK方法的加速效果相對下降。

3 討論

MT正反演中如何選取頻點(diǎn)一直是一個開放的問題[56]，一般采用對數(shù)等間隔選取方法，每個數(shù)量級范圍內(nèi)選取3～4個頻率[57]或者更少[58]，少數(shù)情況下選取5～6個頻率[59]甚至8個頻率[60]。由于采用直接求解方法，每個頻率都需要求解2次大型線性方程組，因此應(yīng)選取盡量少的頻率[56]。

采用本文BRK子空間方法則可以很好地避免上述問題。對于BRK子空間方法，增加求解頻率數(shù)量，只是增加了求解模型降階解方程(式(10))的次數(shù)。由于式(10)是小維度矩陣的傳遞函數(shù)求解，計(jì)算速度很快，因此相應(yīng)增加的計(jì)算時間可以忽略。如果按照一維均勻半空間模型的頻率分布密度采樣，即在1×10-4～1×10Hz頻率范圍內(nèi)計(jì)算51個對數(shù)等間隔分布的頻率，對這些頻點(diǎn)計(jì)算MT響應(yīng)，BRK方法耗時僅為2988.5s，相比21個頻率的MT響應(yīng)計(jì)算時間2894.2s，增加了30個計(jì)算頻率，但求解時間僅增加94.3s，即增加1個頻率所需計(jì)算時間僅為3.1s。而采用直接求解法計(jì)算51個頻率的正演響應(yīng)，所需時間高達(dá)83504.2s，相比計(jì)算21個頻率的計(jì)算時間34391.1s，意味著增加30個頻率所需求解時間增加49113.1s，即增加1個頻率的正演響應(yīng)所需要的時間高達(dá)1637.1s，為BRK方法所需時間的528倍。計(jì)算51個頻率的正演響應(yīng)，BRK方法相比直接求解法的加速比高達(dá)27.9。顯然，求解頻點(diǎn)數(shù)越多，BRK方法的計(jì)算速度優(yōu)勢越明顯。

圖9為采用BRK方法計(jì)算的51個頻率和21個頻率的MT響應(yīng)結(jié)果，這里僅給出通?；趚y和yx分量計(jì)算的視電阻率和相位。如圖9所示，51個頻率的MT視電阻率和相位曲線，相比21個頻率的MT曲線，更能夠體現(xiàn)曲線變化的細(xì)節(jié)特征。因此，在計(jì)算時間允許的條件下，建議選取盡可能多的頻點(diǎn)數(shù)，可得到更精細(xì)的響應(yīng)曲線特征。

圖9 DTM1模型采用BRK方法計(jì)算的不同頻點(diǎn)的視電阻率(上)和相位(下)

4 結(jié)論

本文基于擬態(tài)有限體積離散和塊狀有理Krylov方法實(shí)現(xiàn)了一種求解MT三維正演響應(yīng)的快速算法。通過將源顯式表示為高空中的平面電流源，可將任意頻率的MT正演響應(yīng)表示為一個傳遞函數(shù)與源矢量的乘積。采用塊狀有理Krylov方法，只需要構(gòu)建一次塊狀有理Krylov子空間即可實(shí)現(xiàn)TE和TM極化源在給定頻率范圍內(nèi)所有頻點(diǎn)的MT正演響應(yīng)。采用漸近收斂公式得到塊狀有理Krylov方法的最優(yōu)化單個重復(fù)極點(diǎn)，利用替代算法可快速得到最優(yōu)化的Krylov子空間維度，結(jié)合直接求解器，將MT三維正演的計(jì)算量減小為一次系數(shù)矩陣分解和多次矩陣回代。相對常規(guī)的直接求解方法，本文方法極大地提高了計(jì)算速度。

均勻半空間模型的數(shù)值實(shí)例結(jié)果表明，將平面電流源設(shè)置在空間離散網(wǎng)格中空氣層的最外邊界區(qū)域，能夠最大限度地近似為平面波源，降低正演近似誤差。計(jì)算31個頻率的半空間模型正演響應(yīng)，相比直接求解方法，塊狀有理Krylov方法在保證計(jì)算結(jié)果精度的同時，實(shí)現(xiàn)了22.5倍的加速； 計(jì)算21個頻率的三維DTM1模型的正演響應(yīng)，相比直接求解方法，塊狀有理Krylov方法實(shí)現(xiàn)了11.8倍的加速。求解頻率數(shù)越多，塊狀有理Krylov方法的計(jì)算速度優(yōu)勢越明顯。

本文塊狀有理Krylov方法的主要計(jì)算量為1次系數(shù)矩陣分解和mopt次回代，因此對于中小尺度模型，當(dāng)求解頻率數(shù)量較多時，本文算法在計(jì)算速度上優(yōu)勢顯著。當(dāng)模型的空間離散網(wǎng)格數(shù)量較大、需要求解的頻率數(shù)量較少時，直接分解算法的計(jì)算時間和占用內(nèi)存都會顯著增加，此時基于迭代算法逐次求解各個頻率響應(yīng)的傳統(tǒng)方法在計(jì)算速度上更具有優(yōu)勢。因此，實(shí)際工作中要根據(jù)模型規(guī)模和頻率數(shù)量選擇合適的正演求解方法。