周建美 劉文韜 魯凱亮 李 貅

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

0 引言

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

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

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

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

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

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

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

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

1 正演算法

1.1 控制方程離散

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

?×E+iωB=0

(1)

?×μ-1B-σE=S

(2)

式中：E和B分別表示電場(chǎng)強(qiáng)度矢量和磁場(chǎng)強(qiáng)度矢量； 角頻率ω=2πf，其中f是頻率；μ是磁導(dǎo)率；σ是電導(dǎo)率；S是外加源項(xiàng)，顯式表示為平面電流分布[46]。采用自然邊界條件，利用基于交錯(cuò)網(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分別為離散形式的電場(chǎng)、磁場(chǎng)和發(fā)射源項(xiàng)。對(duì)于三維MT正演，同時(shí)考慮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)中的磁場(chǎng)，得到關(guān)于電場(chǎng)的控制方程

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

(8)

e(ω)=gω(A)X

(9)

1.2 塊狀有理Krylov方法

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

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

(10)

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

(11)

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

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

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

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

1.3 MT響應(yīng)

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

(12)

其中下標(biāo)“1”和“2”分別代表TE和TM極化模式。基于阻抗張量Z定義視電阻率和相位[13]分別為

(13)

(14)

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

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

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

2.1 均勻半空間模型

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

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

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

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

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

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

為了分析BRK方法的計(jì)算效率，對(duì)比本文BRK方法、非塊狀有理Krylov方法[42](RK)和逐次采用直接求解器PARDISO[53]求解各個(gè)頻率的正演響應(yīng)的直接求解法(Direct)的計(jì)算結(jié)果。雖然一維模型的TE模式和TM模式的響應(yīng)是相同的，但是為了通用性考慮，在比較算法計(jì)算效率時(shí)，本文依然以計(jì)算TE和TM模式響應(yīng)的總計(jì)算時(shí)間為依據(jù)。上述三種算法計(jì)算的MT響應(yīng)基本相同，因此本文僅對(duì)計(jì)算時(shí)間進(jìn)行比較，結(jié)果如表3所示?？梢?，本文BRK方法只需1次系數(shù)矩陣分解及70次矩陣回代； RK方法分別求解TE和TM模式的響應(yīng)，因此需要1次系數(shù)矩陣分解和140次矩陣回代； 直接求解法(Direct)需對(duì)每個(gè)頻率進(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ì)算量及耗時(shí)對(duì)比

2.2 三維DTM1模型

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

DTM1模型的均勻半空間包含3個(gè)塊狀電性異常體，具體參數(shù)見表4和圖6。坐標(biāo)原點(diǎn)位于塊體1中心點(diǎn)的正上方地表，觀測(cè)點(diǎn)位于坐標(biāo)原點(diǎn)，fmin=10-4Hz，fmax=10Hz，在計(jì)算頻率范圍內(nèi)對(duì)數(shù)等間隔取21個(gè)頻點(diǎn)。計(jì)算區(qū)域離散為72×74×72個(gè)網(wǎng)格，最小網(wǎng)格長(zhǎ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子空間解的相對(duì)誤差error隨m的變化曲線。相比一維模型，三維DTM1模型的求解頻率范圍更大，因此最小頻率fmin對(duì)應(yīng)的誤差收斂率更小，從而相對(duì)誤差穩(wěn)定時(shí)所對(duì)應(yīng)的m更大。由圖7可知，m>150時(shí)，相對(duì)誤差基本恒定，因此選取mopt=150。

圖7 fmin=1×10-4Hz時(shí)有理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)由于視電阻率值太小沒有顯示，相位分布也相對(duì)雜亂； 但是在低頻段(<1Hz)，本文計(jì)算結(jié)果與積分方程法和有限元法的結(jié)果基本一致。因此，本文BRK方法能夠較精確地正演MT三維復(fù)雜模型的電磁響應(yīng)。

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

3 討論

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

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

圖9為采用BRK方法計(jì)算的51個(gè)頻率和21個(gè)頻率的MT響應(yīng)結(jié)果，這里僅給出通?；趚y和yx分量計(jì)算的視電阻率和相位。如圖9所示，51個(gè)頻率的MT視電阻率和相位曲線，相比21個(gè)頻率的MT曲線，更能夠體現(xiàn)曲線變化的細(xì)節(jié)特征。因此，在計(jì)算時(shí)間允許的條件下，建議選取盡可能多的頻點(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)表示為一個(gè)傳遞函數(shù)與源矢量的乘積。采用塊狀有理Krylov方法，只需要構(gòu)建一次塊狀有理Krylov子空間即可實(shí)現(xiàn)TE和TM極化源在給定頻率范圍內(nèi)所有頻點(diǎn)的MT正演響應(yīng)。采用漸近收斂公式得到塊狀有理Krylov方法的最優(yōu)化單個(gè)重復(fù)極點(diǎn)，利用替代算法可快速得到最優(yōu)化的Krylov子空間維度，結(jié)合直接求解器，將MT三維正演的計(jì)算量減小為一次系數(shù)矩陣分解和多次矩陣回代。相對(duì)常規(guī)的直接求解方法，本文方法極大地提高了計(jì)算速度。

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

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