劉寄仁， 湯井田,2,3,4*， 任政勇,2,3,4， 張繼鋒

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長(zhǎng)沙 410083 2 中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)沙 410083 3 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)沙 410083 4 自然資源部覆蓋區(qū)深部資源勘查工程技術(shù)創(chuàng)新中心， 合肥 230001 5 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院， 西安 710061

0 引言

陸地可控源電磁目前被廣泛應(yīng)用于礦產(chǎn)資源、天然氣、煤田采空區(qū)富水性調(diào)查等環(huán)境與工程地球物理領(lǐng)域(湯井田和何繼善, 2005).正演模擬可以為研究實(shí)際地質(zhì)響應(yīng)特征提供參照，因此需要開發(fā)快速的三維正演計(jì)算方法.

陸地可控源電磁的主流正演方法有積分方程法(Avdeev et al., 2002; Zhdanov et al., 2006; 任政勇等, 2017; 湯井田等, 2018)、有限差分法(Mackie et al., 1994; Newman and Alumbaugh, 2002; Streich, 2009; Malovichko et al., 2019)、有限體積法(Jahandari and Farquharson, 2014, 2015; Peng et al., 2018; Liu et al., 2019)、有限單元法(Mitsuhata and Uchida, 2004; Schwarzbach et al., 2011; Ren et al., 2013, 2014; Grayver and Kolev, 2015; Um et al., 2017; Liu et al., 2018).有限單元法近年來(lái)得到了廣泛的應(yīng)用，主要因?yàn)樗哂欣碚擉w系成熟、網(wǎng)格剖分靈活等特點(diǎn)，結(jié)合非結(jié)構(gòu)化網(wǎng)格可以模擬復(fù)雜地形及地質(zhì)結(jié)構(gòu)，便于在場(chǎng)源、測(cè)點(diǎn)及電性分界面處進(jìn)行網(wǎng)格加密，從而提高數(shù)值計(jì)算的精度.有限單元法正演可以采用總場(chǎng)公式和二次場(chǎng)公式.總場(chǎng)公式需要對(duì)場(chǎng)源附近的電磁場(chǎng)做精細(xì)處理，并且在計(jì)算邊界上需要施加精細(xì)的邊界條件，但總場(chǎng)公式能夠有效模擬任意起伏地形(李建慧等, 2016; 邱長(zhǎng)凱等, 2018).二次場(chǎng)公式的模擬變量為散射場(chǎng)，散射場(chǎng)在源附近一般不再具有奇異性，因此具有模擬精度高的特點(diǎn)，但是二次場(chǎng)公式需要背景模型具有解析的一次場(chǎng)表達(dá)式(張林成等, 2017).一般來(lái)說(shuō), 地形背景模型不具備解析的一次場(chǎng)表達(dá)式，二次場(chǎng)公式處理地形存在困難.因此，本文選擇總場(chǎng)公式作為開發(fā)可處理任意起伏地形情況的陸地可控源電磁快速正演方法.

陸地可控源電磁有限元正演最終形成一個(gè)大型的復(fù)系數(shù)方程組.該方程組和頻率有關(guān)，各頻點(diǎn)之間相互獨(dú)立，可以使用基于頻點(diǎn)的CPU并行求解技術(shù)來(lái)實(shí)現(xiàn)加速，但該加速方案受限于計(jì)算機(jī)硬件設(shè)備.另一種方案是對(duì)初始線性方程組的系數(shù)矩陣進(jìn)行降階，在保留系統(tǒng)性質(zhì)的基礎(chǔ)上，通過(guò)降階手段，使用一個(gè)較小的系統(tǒng)來(lái)近似該矩陣，從而得到一個(gè)近似解.Krylov子空間投影算法是一類典型的模型降階算法，通過(guò)構(gòu)建Krylov子空間的正交基，將初始矩陣投影到子空間上，得到一個(gè)較小維度的降階矩陣，使用該降階矩陣來(lái)替代初始矩陣，便可以得到初始線性方程組的近似解.早期，基于m維多項(xiàng)式Krylov子空間的SLDM(spectral Lanczos decomposition method)方法被應(yīng)用到計(jì)算電磁領(lǐng)域(Druskin and Knizhnerman, 1988, 1994)，但是這種多項(xiàng)式Krylov方法不穩(wěn)定，收斂速度依賴于頻率范圍和模型的電導(dǎo)率反差，需要一個(gè)較大的m維子空間才能保證誤差收斂到較小的水平，同時(shí)隨著子空間維度m的增大，誤差會(huì)進(jìn)一步積累.為解決大型稀疏矩陣的特征值問(wèn)題，Ruhe (1984,1994)提出了基于有理函數(shù)的Krylov子空間方法，與多項(xiàng)式Krylov子空間方法相比，該方法更加穩(wěn)定，所需的子空間維度較小，能夠減小計(jì)算量.

在地球物理領(lǐng)域中，Druskin等(2009)和Knizhnerman等(2009)通過(guò)誤差分析理論，分別提出了時(shí)間域和頻率域的有理Krylov子空間的多極點(diǎn)選擇方案，其優(yōu)點(diǎn)是收斂誤差不依賴于模型的電導(dǎo)率反差，但對(duì)于m維的子空間，需要進(jìn)行m次矩陣方程求解，如何高效的求解這些方程直接決定了有理Krylov子空間模型降階算法的效率.為減少極點(diǎn)的個(gè)數(shù)，加快子空間的構(gòu)建，B?rner等(2008)通過(guò)選擇較少的參考極點(diǎn)構(gòu)建有理Krylov子空間，使用間接法求解時(shí)間域瞬變電磁問(wèn)題，在一定精度條件下，大大提高了求解效率.B?rner等(2015)提出了替代最優(yōu)化算法，給出了一定時(shí)間范圍內(nèi)的極點(diǎn)選擇方案，使用基于非結(jié)構(gòu)化網(wǎng)格的矢量有限元方法實(shí)現(xiàn)了時(shí)間域瞬變電磁的快速正演求解.Qiu等(2019)對(duì)于寬時(shí)間范圍的雙極點(diǎn)選擇方案進(jìn)行了修正，使用塊Krylov方法，實(shí)現(xiàn)了時(shí)間域海洋可控源電磁的快速正演計(jì)算，有效的提高了多源正演模擬的速度.為實(shí)現(xiàn)頻率域可控源電磁法的快速正演求解，周建美等(2018)使用單極點(diǎn)模型降階算法，結(jié)合擬態(tài)有限體積進(jìn)行了三維海洋可控源電磁法的模擬.張繼鋒等(2020)使用單極點(diǎn)模型降階算法，結(jié)合六面體節(jié)點(diǎn)有限元進(jìn)行了陸地CSAMT的正演求解，但是誤差較大，且六面體網(wǎng)格不適合模擬復(fù)雜地形和不規(guī)則地質(zhì)體.此外，正交化算法也是影響計(jì)算效率的重要因素之一.在構(gòu)建Krylov子空間的正交基時(shí)，常使用Arnoldi正交化算法(B?rner et al., 2008, 2015;周建美等, 2018; 張繼鋒等, 2020)，雖然該算法較穩(wěn)定，但是當(dāng)子空間的維數(shù)增加時(shí)，需要進(jìn)行大量的大型矩陣向量乘積運(yùn)算，乘積運(yùn)算的次數(shù)隨著子空間的維度數(shù)呈平方增加，嚴(yán)重影響計(jì)算效率.而Lanczos算法的矩陣向量乘積的次數(shù)隨著子空間維度呈線性增長(zhǎng)，因此使用Lanczos算法能夠更加快速的構(gòu)建Krylov子空間的正交基.

為加快多頻正演的求解速度，本文結(jié)合基于非結(jié)構(gòu)化網(wǎng)格的矢量有限元方法和單極點(diǎn)模型降階算法，使用Lanczos算法構(gòu)建子空間的正交基，實(shí)現(xiàn)了三維陸地多頻可控源電磁的快速正演.首先，采用非結(jié)構(gòu)化網(wǎng)格進(jìn)行空間離散，得到了有限元線性方程組.然后選擇一定頻率范圍內(nèi)的單個(gè)最優(yōu)化極點(diǎn)，將大型稀疏復(fù)系數(shù)方程組的求解問(wèn)題，轉(zhuǎn)化為實(shí)數(shù)方程組求解，使用Lanczos算法高效的構(gòu)建了有理Krylov子空間的正交基，求得初始矩陣在Krylov子空間上的正交投影.投影矩陣能適應(yīng)所有頻率的計(jì)算，其維數(shù)是遠(yuǎn)小于初始矩陣的，因此使用投影矩陣來(lái)替代初始矩陣，可以實(shí)現(xiàn)矩陣方程的快速求解.最后通過(guò)一維層狀模型和三維塊狀體模型進(jìn)行了加速比測(cè)試，給出了本文算法的精度對(duì)比，并且詳細(xì)的給出了整個(gè)過(guò)程的運(yùn)算時(shí)間對(duì)比和內(nèi)存消耗對(duì)比，驗(yàn)證了本文算法具有計(jì)算資源占用率低的優(yōu)勢(shì).

1 方法

1.1 CSEM有限元線性方程組

設(shè)時(shí)間諧變因子為eiω t，電場(chǎng)強(qiáng)度E的微分控制方程為：

(1)

為求解式(1)，B?rner等(2008)使用基于非結(jié)構(gòu)化網(wǎng)格的矢量有限元法結(jié)合有理Krylov子空間模型降階算法得到了近似解，但該方法在高頻和低頻部分誤差較大，這雖然可以通過(guò)調(diào)整極點(diǎn)的位置和個(gè)數(shù)來(lái)改善，但無(wú)疑會(huì)增加矩陣分解的次數(shù)，從而增加計(jì)算消耗.周建美等(2018)和張繼鋒等(2020)通過(guò)單極點(diǎn)模型降階算法，分別使用擬態(tài)有限體積法和六面體節(jié)點(diǎn)有限元法實(shí)現(xiàn)了一定頻率范圍內(nèi)的快速求解.本文使用基于非結(jié)構(gòu)化網(wǎng)格的矢量有限元法和單極點(diǎn)模型降階算法來(lái)對(duì)式(1)進(jìn)行快速求解.

首先，采用開源的非結(jié)構(gòu)化四面體網(wǎng)格剖分器Tetgen(Si, 2015)將計(jì)算區(qū)域剖分成多個(gè)互斥的四面體單元.然后,以圖1所示的四面體單元為例，單元內(nèi)部的電場(chǎng)值可以通過(guò)如下的線性插值公式表示:

圖1 四面體單元Fig.1 Tetrahedral element

(2)

(3)

(4)

將(2)式代入(4)式，利用矢量恒等變換，并假設(shè)網(wǎng)格邊界取得足夠遠(yuǎn)，忽略邊界積分項(xiàng)，可得：

(5)

對(duì)于任意單元e，可將式(5)寫成矩陣的形式：

Re=(Ce+iωMe)Ee+iωbe,

(6)

(C+iωM)E=-iωb,

(7)

(7)式便是有限元分析最終得到的大型線性方程組，其系數(shù)矩陣為大型復(fù)數(shù)稀疏矩陣.對(duì)于多頻點(diǎn)問(wèn)題，必須進(jìn)行多次矩陣方程的求解，需要消耗大量的計(jì)算內(nèi)存和運(yùn)算時(shí)間，嚴(yán)重影響正演計(jì)算的效率.

1.2 基于模型降階的多頻線性方程組加速

為了提高多頻可控源電磁的正演效率，本文引入一種不依賴于頻點(diǎn)個(gè)數(shù)的有理Krylov子空間模型降階算法，在一定精度條件下，該算法只需要一次矩陣分解，便可對(duì)系統(tǒng)矩陣實(shí)現(xiàn)降階，從而達(dá)到對(duì)方程組進(jìn)行快速求解的目的(Güttel, 2010, 2013; Zhou et al., 2018b, 2020;周建美等, 2018; 張繼鋒等, 2020).對(duì)式(7)做一個(gè)簡(jiǎn)單的變形，令A(yù)=M-1C，X=M-1b，方程的解E可寫為：

E=-iω(A+iωI)-1X=f(A)X.

(8)

為解決此類形如f(A)X的問(wèn)題，構(gòu)建有理Krylov子空間為：

(9)

(10)

(11)

E≈‖X‖MVm+1f(Am+1)e1,

(12)

其中e1=[1,0,0,…,0]T∈m+1.可以發(fā)現(xiàn)，一旦求得降階矩陣Am+1，那么f(Am+1)是容易求解的，對(duì)于多頻問(wèn)題，求解式(12)是快速的.

構(gòu)建上述有理Krylov子空間的正交基的方法有兩種，一種是Arnoldi正交化算法，該算法較穩(wěn)定，但是每生成一個(gè)正交基向量需要和已生成的所有基向量進(jìn)行正交，我們?cè)谑?10)中已經(jīng)給出了構(gòu)建正交基的遞歸表達(dá)式.另一種方法是Lanczos算法，雖然Lanczos算法可能會(huì)不穩(wěn)定，但仍然會(huì)對(duì)初始矩陣形成有效的近似(Druskin and Remis, 2013)，每生成一個(gè)正交基向量只需利用已生成的兩個(gè)基向量進(jìn)行運(yùn)算.其正交基的遞歸表達(dá)式為：

vj+1=((A-ξjI)-1vj-αjvj-βj-1vj-1)/βj

=((C-ξjM)-1Mvj-αjvj-βj-1vj-1)/βj,(13)

其中αj=(Mvj)T(C-ξjM)-1Mvj，β0v0=0，βj=‖(C-ξjM)-1Mvj-αjvj-βj-1vj-1‖M.利用M正交投影，同樣可以得到式(8)的近似表達(dá)式(12).相對(duì)Arnoldi算法而言，Lanczos算法中大型矩陣向量的乘積運(yùn)算次數(shù)大大減少，可以提高構(gòu)建正交基的效率，因此本文選擇Lanczos算法來(lái)構(gòu)建正交基，算法的詳細(xì)過(guò)程在表1中給出.

表1 基于M正交的Lanczos算法Table 1 Lanczos algorithm based on M-orthogonal

從表1可以看出，模型降階算法將原來(lái)的復(fù)系數(shù)方程組求解轉(zhuǎn)換為實(shí)系數(shù)方程組的求解，這顯然可以提高運(yùn)算效率和節(jié)省內(nèi)存，但正交基的構(gòu)建需要選擇不同的實(shí)數(shù)極點(diǎn)ξj，這帶來(lái)的一個(gè)問(wèn)題便是需要進(jìn)行多次矩陣方程的求解.位移求逆技術(shù)(Moret and Novati, 2004)是一種特殊的有理Krylov子空間模型降階方法，其僅需選擇一個(gè)重復(fù)極點(diǎn)滿足：ξ0=ξj<0,ξ0?Λ(A).這種單極點(diǎn)模型降階算法只需要一次Cholesky分解，多次矩陣回代，便可快速構(gòu)建子空間.因此本文選擇單極點(diǎn)模型降階算法來(lái)減少矩陣方程的求解次數(shù)，只需在表1中第03步j(luò)=1時(shí)做一次Cholesky分解，并將矩陣分解的結(jié)果保存，之后的循環(huán)只需調(diào)用即可.

(14)

2 數(shù)值實(shí)驗(yàn)

本文進(jìn)行數(shù)值模擬的計(jì)算環(huán)境為：Ubuntu 18.04和Intel Visual Fortran 19.1,硬件為：Intel(R) Core(TM) i7-9750H CPU @2.60 GHz 2.59 GHz，16 GB個(gè)人PC機(jī).

2.1 加速比測(cè)試

2.1.1 一維模型

為驗(yàn)證本文算法的正確性和高效性，首先構(gòu)建一維層狀模型，第一層電阻率為100 Ωm，厚度為1000 m，第二層電阻率為500 Ωm.場(chǎng)源的中心設(shè)置在坐標(biāo)原點(diǎn)，沿著x軸放置，長(zhǎng)度為100 m，電流大小為10 A.發(fā)射頻率為1～10000 Hz對(duì)數(shù)等間隔分布的100個(gè)頻點(diǎn).使用Tetgen(Si, 2015)進(jìn)行網(wǎng)格剖分，未知數(shù)為629370.

使用本文算法進(jìn)行正演計(jì)算(3DMOR_Lanczos)，同一維解析解(1D solution)、常規(guī)有限元求解(3DFEM)、基于Arnoldi方法的模型降階求解(3DMOR_Arnoldi)進(jìn)行精度和時(shí)間對(duì)比.3DMOR_Lanczos和3DMOR_Arnoldi構(gòu)建Krylov子空間的維數(shù)均為120.圖2給出了測(cè)點(diǎn)為(0 m,-4000 m,0 m)處電場(chǎng)Ex的響應(yīng)曲線及相對(duì)誤差曲線.從圖中可以看出，3DMOR_Lanczos同其他兩種方法的Ex響應(yīng)曲線均吻合較好，同一維解析解之間的最大誤差不超過(guò)3%，驗(yàn)證了本文算法的正確性.3DMOR_Lanczos同3DFEM、3DMOR_Arnoldi之間的相對(duì)誤差很小，說(shuō)明3DMOR _Lanczos對(duì)3DFEM的近似效果較好，且與Arnoldi方法的結(jié)果吻合.在y軸上布置一條測(cè)線，起始位置為(0 m,-5000 m,0 m)，終止位置為(0 m,-100 m,0 m)，圖3中呈現(xiàn)的是105 Hz時(shí)測(cè)線上的Ex響應(yīng)的變化.從中可知，3DMOR_Lanczos同1D solution之間的最大誤差小于3%，同3DFEM、3DMOR_Arnoldi之間的誤差很小，進(jìn)一步驗(yàn)證了本文算法的正確性.表2給出了計(jì)算資源消耗的對(duì)比.加速比定義為其他方法與本文算法(3DMOR_Lanczos)的運(yùn)算時(shí)間的比值.本文算法將復(fù)數(shù)矩陣轉(zhuǎn)變?yōu)閷?shí)數(shù)矩陣求解，因此與3DFEM相比，3DMOR_Lanczos消耗的內(nèi)存更少，且計(jì)算時(shí)間大幅減少，加速比達(dá)到了24.9.與3DMOR_Arnoldi相比，二者消耗的內(nèi)存相近，但加速比達(dá)到了2.7，表明了本文算法的高效率.

表2 運(yùn)算時(shí)間對(duì)比Table 2 Comparison of calculation time

表3給出了3DMOR_Lanczos和3DMOR_Arnoldi的運(yùn)行時(shí)間分布，由于矩陣的規(guī)模較小，結(jié)合Lanczos算法的高效性，Pardiso矩陣分解和構(gòu)建正交基消耗的時(shí)間占比僅為40.6%，最終的頻點(diǎn)運(yùn)算消耗了大部分的運(yùn)算時(shí)間.3DMOR_Lanczos和3DMOR_Arnoldi兩種算法的唯一不同在于正交基的構(gòu)建方法不一樣，我們注意到，3DMOR_Arnoldi的正交基構(gòu)建占據(jù)了整個(gè)過(guò)程的66.39%，這說(shuō)明使用Lanczos算法構(gòu)建正交基更加高效.為方便說(shuō)明和閱讀，后文將3DMOR_Lanczos縮寫為3DMOR.

表3 3DMOR的運(yùn)算時(shí)間分布Table 3 Calculation time distribution of 3DMOR

圖2 3DMOR_Lanczos同1D solution、3DFEM、3DMOR_Arnoldi的測(cè)深曲線對(duì)比(a) 不同算法的Ex響應(yīng)對(duì)比圖; (b) 3DMOR_Lanczos同其他方法之間的相對(duì)誤差.Fig.2 Comparison of sounding curves of 3DMOR_Lanczos with 1D solution, 3DFEM and 3DMOR_Arnoldi(a) Comparison of Ex responses of different algorithms; (b) The relative error between 3DMOR_Lanczos and other methods.

圖3 頻率為105 Hz時(shí), 3DMOR_Lanczos同1D solution、3DFEM、3DMOR_Arnoldi對(duì)比(a) 不同算法的Ex響應(yīng)對(duì)比圖; (b) 3DMOR_Lanczos同其他方法之間的相對(duì)誤差.Fig.3 Comparison of 3DMOR_Lanczos with 1D solution, 3DFEM and 3DMOR_Arnoldi at 105 Hz(a) Comparison of Ex responses of different algorithms; (b) The relative error between 3DMOR_Lancozs and other methods.

2.1.2 三維模型

為進(jìn)一步說(shuō)明本文算法的高效性和對(duì)三維模型的適應(yīng)性，設(shè)計(jì)如圖4所示的三維模型進(jìn)行對(duì)比測(cè)試.圖4a為模型的剖面圖，圖4b為模型的平面圖，圖中虛線為測(cè)線.在電阻率為50 Ωm的均勻半空間中嵌入一個(gè)大小為120 m×200 m×400 m的低阻異常體，其中心坐標(biāo)為(1000 m,0 m,300 m)，電阻率為5 Ωm.100 m長(zhǎng)的源沿著x軸放置，源的中心坐標(biāo)為(50 m,0 m,0 m)，發(fā)射電流為1 A，發(fā)射頻率為1～10000 Hz對(duì)數(shù)等間隔分布的32個(gè)頻點(diǎn).網(wǎng)格離散所形成的未知數(shù)為889023.

圖4 三維低阻體模型示意圖(a) 模型的剖面圖; (b) 模型的平面圖.Fig.4 Sketch of 3D low resistance body(a) Section view of the model; (b) Plan view of the model.

將本文數(shù)值解(3DMOR)與常規(guī)有限元方法(3DFEM)、基于磁矢勢(shì)的有限元方法(FE)(Ansari and Farquharson, 2014)、積分方程法(IE)(Farquharson and Oldenburg, 2002; Zhou et al., 2018a)的數(shù)值解進(jìn)行了Ex實(shí)部和虛部響應(yīng)曲線的對(duì)比，正演計(jì)算頻率為3 Hz，如圖5a、圖6a所示.從圖中可知，由于地下低阻體的存在，地下電流被低阻體“吸引”，導(dǎo)致地表的電流密度降低，從而出現(xiàn)了電場(chǎng)Ex響應(yīng)的實(shí)部和虛部曲線在x為900～1100 m范圍內(nèi)均呈現(xiàn)向下凹陷現(xiàn)象.圖5b和圖6b給出了Ex響應(yīng)的相對(duì)誤差曲線，除IE方法外，本文算法同其他兩種方法之間的相對(duì)誤差未超過(guò)3%，這說(shuō)明利用本文算法進(jìn)行三維模型的正演是可行的.

圖5 頻率為3 Hz時(shí), 3DMOR同多種數(shù)值模擬方法的電場(chǎng)實(shí)部的對(duì)比(a) 不同算法的Ex實(shí)分量對(duì)比圖; (b) 3DMOR同其他方法之間的相對(duì)誤差.Fig.5 Comparison of the real part of electric field between 3DMOR and other numerical methods at 3 Hz(a) Comparison of the real part of Ex responses of different algorithms; (b) The relative error between 3DMOR and other methods.

圖6 頻率為3 Hz時(shí), 3DMOR同多種數(shù)值模擬方法的電場(chǎng)虛部的對(duì)比(a) 不同算法的Ex虛分量對(duì)比圖; (b) 3DMOR同其他方法之間的相對(duì)誤差.Fig.6 Comparison of the imaginary part of electric field between 3DMOR and other numerical methods at 3 Hz(a) Comparison of the imaginary part of Ex responses of different algorithms; (b) The relative error between 3DMOR and other methods.

圖7給出了測(cè)點(diǎn)(1000 m,0 m,0 m)處3DMOR和3DFEM的Ex響應(yīng)曲線及相對(duì)誤差曲線，可以看出3DMOR對(duì)3DFEM的近似效果較好.同樣在表4給出了計(jì)算資源消耗的對(duì)比，從中可知加速比達(dá)到了17.6.與表2相比可知，隨著未知數(shù)的增加，3DMOR更節(jié)省內(nèi)存.通過(guò)表3我們可以發(fā)現(xiàn)，隨著未知數(shù)的增加，矩陣分解和正交化過(guò)程將消耗大部分的運(yùn)算時(shí)間，因此本文采用的高效的Lanczos正交化算法是可以提升正演計(jì)算效率的.

圖7 中心測(cè)點(diǎn)的測(cè)深曲線，3DMOR同3DFEM的對(duì)比(a) 兩種算法的Ex響應(yīng)對(duì)比圖; (b) 3DMOR同3DFEM之間的相對(duì)誤差.Fig.7 Sounding curve of central measuring point, comparison between 3DMOR and 3DFEM(a) Comparison of Ex responses of two algorithms; (b) The relative error between 3DMOR and 3DFEM.

表4 運(yùn)算時(shí)間對(duì)比Table 4 Comparison of calculation time

2.2 頻率分段策略

為適應(yīng)更寬頻的可控源電磁的正演計(jì)算，首先需要探究頻帶寬度和Krylov子空間維度對(duì)求解精度和速度的影響.基于上文的1D層狀模型，固定最高頻率為10000 Hz，保持正演計(jì)算頻率的數(shù)量不變，逐步改變頻率的范圍和Krylov子空間維度m，fmax表示最高頻率，fmin表示最低頻率，利用式(14)得到了圖8所示的誤差曲線和運(yùn)算時(shí)間分布圖.從圖8a可以看出，當(dāng)頻帶范圍越小時(shí)，整體誤差收斂的速度越快，只需要構(gòu)建較小的Krylov子空間便能達(dá)到較高的精度.而從圖8b所示的計(jì)算時(shí)間分布上也可以看出，由于只改變了Krylov子空間的維度，因此四條曲線的變化趨勢(shì)是吻合的.隨著Krylov子空間維度的增加，運(yùn)行時(shí)間基本呈現(xiàn)線性增長(zhǎng)，這是Lanczos算法的計(jì)算優(yōu)勢(shì).

圖8 不同頻帶寬度下，3DMOR的(a)整體誤差和(b)運(yùn)行時(shí)間隨Krylov子空間維度的變化Fig.8 Under different bandwidth, (a) the overall error and (b) the running time of 3DMOR change with the dimension of Krylov subspace

在頻率域可控源電磁的正演計(jì)算中，當(dāng)發(fā)射頻率的范圍較大時(shí)，為控制高頻和低頻的精度，不僅要求淺部的網(wǎng)格剖分更精細(xì)，還要求網(wǎng)格邊界取得足夠遠(yuǎn).特別是對(duì)于層狀介質(zhì)而言，網(wǎng)格單元數(shù)量增長(zhǎng)很快，會(huì)浪費(fèi)一定的計(jì)算資源.在使用3DMOR進(jìn)行寬頻正演計(jì)算時(shí)，基于Krylov子空間維度較大和網(wǎng)格單元較多這兩個(gè)特征，本文提出頻率分段策略：針對(duì)寬頻問(wèn)題(log(fmax/fmin)>4)，對(duì)頻率進(jìn)行分段，同時(shí)輸入兩套網(wǎng)格，使用并行算法進(jìn)行3DMOR的運(yùn)算.這不僅可以減小矩陣方程的規(guī)模，還可以減小Krylov子空間的維度，加快誤差的收斂速度.

對(duì)于上述的1D層狀模型，擴(kuò)大頻率計(jì)算范圍為0.01～10000 Hz，正演計(jì)算100個(gè)頻率，測(cè)點(diǎn)坐標(biāo)為(0 m,-4000 m,0 m).改變Krylov子空間的維數(shù)m，其余參數(shù)保持不變.圖9給出了四種不同的Krylov子空間維度下的Ex響應(yīng)曲線.可以發(fā)現(xiàn)3DMOR在高頻段的Ex響應(yīng)出現(xiàn)了波動(dòng),m需要達(dá)到200才能和1D solution 曲線呈現(xiàn)較好的吻合，此時(shí)的運(yùn)算時(shí)間為240 s.雖然隨著Krylov子空間維度的增加，誤差會(huì)逐步減小，但無(wú)論是選擇式(10)所示的Arnoldi算法還是式(13)描述的Lanczos算法，過(guò)大的子空間需求都會(huì)增加構(gòu)建正交基的計(jì)算量，影響計(jì)算效率.

圖9 當(dāng)頻率計(jì)算范圍為[0.01,10000]Hz時(shí)，3DMOR采取不同Krylov子空間維度的測(cè)深曲線同1D solution對(duì)比Fig.9 When the frequency calculation range is [0.01,10000] Hz, the sounding curves of 3DMOR with different Krylov subspace dimensions are compared with 1D solution

我們使用頻率分段策略來(lái)解決這一問(wèn)題.對(duì)于高頻部分，設(shè)置網(wǎng)格邊界為±20 km，方程未知數(shù)為356637，構(gòu)建Krylov子空間的維度為80.對(duì)于低頻部分，設(shè)置網(wǎng)格邊界為±150 km，方程未知數(shù)為348119，同樣構(gòu)建Krylov子空間的維度為80.最終總內(nèi)存消耗為3.33 GB，總運(yùn)行時(shí)間為101 s.圖10給出了進(jìn)行頻率分段之后的運(yùn)行結(jié)果，黑色的虛線為頻率分段的交界點(diǎn)，結(jié)果顯示最大誤差未超過(guò)3%，表明了本文提出的頻率分段策略對(duì)于寬頻可控源電磁的正演計(jì)算是可行的.

圖10 當(dāng)頻率計(jì)算范圍為[0.01,10000]Hz時(shí)，采取頻率分段策略的計(jì)算結(jié)果(a) 3DMOR同1D solution的測(cè)深曲線對(duì)比; (b) 相對(duì)誤差.Fig.10 Calculation results of frequency segmentation strategy when frequency calculation range is [0.01,10000]Hz(a) Comparison of sounding curve between 3DMOR and 1D solution; (b) Relative error.

2.3 三維地質(zhì)體響應(yīng)特征分析

設(shè)計(jì)如圖11所示的三維低阻球體模型，圖12為網(wǎng)格剖分示意圖.球體的半徑為200 m，中心埋深位于(0 m,5000 m,400 m)，電阻率為10 Ωm，圍巖電阻率為100 Ωm.x方向的電流源放置在坐標(biāo)原點(diǎn)，其長(zhǎng)度為200 m，發(fā)射電流為10 A，發(fā)射頻率為0.01～10000 Hz對(duì)數(shù)等間隔分布的100個(gè)頻點(diǎn).本例在地表布置了11條測(cè)線，線距為200 m，y坐標(biāo)范圍為-1000～1000 m.每條測(cè)線長(zhǎng)度為2000 m，含11個(gè)測(cè)點(diǎn)，點(diǎn)距為200 m，x坐標(biāo)范圍為4000～6000 m.仍然采用上文的頻率分段策略來(lái)進(jìn)行正演求解.對(duì)于高頻段，設(shè)置網(wǎng)格邊界為±20 km，方程未知數(shù)為578603，構(gòu)建Krylov子空間的維度為80.對(duì)于低頻段，網(wǎng)格邊界為±150 km，方程未知數(shù)為387979，構(gòu)建Krylov子空間的維度為80.兩套網(wǎng)格并行計(jì)算，最終消耗總內(nèi)存為4.11 GB，總運(yùn)行時(shí)間為148 s.

圖11 三維低阻球體模型示意圖(a) 模型的平面圖； (b) 模型的切片圖.Fig.11 Sketch of 3D low resistance sphere model(a) Plan view of the model; (b) Section view of the model.

圖12 網(wǎng)格剖分示意圖Fig.12 Schematic diagram of mesh generation

圖13給出了四個(gè)頻率的磁場(chǎng)Hy分量、電場(chǎng)Ex分量以及廣域視電阻率的切片圖，四個(gè)頻率分別為10000 Hz、100 Hz、1 Hz、0.01 Hz，其中廣域視電阻率利用Ex求得(李帝銓等, 2013; He, 2018;武建平等, 2020；Yu et al., 2020).從圖13(a,b,c,d)中可以看出，磁場(chǎng)Hy分量對(duì)異常體的響應(yīng)較微弱，無(wú)法從切片圖中觀察到異常場(chǎng)造成的畸變.如圖13(f,g,h,j,k,l)所示，由于低阻異常體的存在，使得地表電流密度減小，因此電場(chǎng)Ex分量和廣域視電阻率均出現(xiàn)凹陷，在頻率為100 Hz、1 Hz、0.01 Hz均有不同程度的響應(yīng)，能夠較好的反映異常體的位置和屬性.圖14給出的是兩條中心測(cè)線(x=0 m測(cè)線,y=5000 m測(cè)線)的廣域視電阻率擬斷面圖，可以看出，在高頻部分廣域視電阻率趨向于背景電阻率，在低頻部分呈現(xiàn)低阻響應(yīng)形態(tài).

圖13 頻率分別為10000 Hz、100 Hz、1 Hz、0.01 Hz的Hy、Ex、廣域視電阻率切片圖(a)—(d)磁場(chǎng)Hy分量切片圖; (e)—(h)電場(chǎng)Ex分量切片圖; (i)—(l)廣域視電阻率切片圖.Fig.13 Slices of Hy, Ex and wide field apparent resistivity with frequencies of 10000 Hz, 100 Hz, 1 Hz and 0.01 Hz respectively(a)—(d) Slices of magnetic field Hy; (e)—(h) Slices of electric field Ex; (i)—(l) Slices of wide field apparent resistivity.

圖14 x=0測(cè)線和y=5000測(cè)線的廣域視電阻率擬斷面圖Fig.14 Pseudo-section of wide field apparent resistivity with line x=0 and line y=5000

圖15給出了三個(gè)測(cè)點(diǎn)處的測(cè)深曲線，測(cè)點(diǎn)分別位于(-1000 m,5000 m,0 m)、(1000 m,5000 m,0 m)、(0 m,5000 m,0 m)，圖15a為電場(chǎng)Ex分量的變化曲線，圖15b為廣域視電阻率變化曲線，黑色的虛線為頻率分段的交界點(diǎn).從圖中可以看出，隨著頻率的減小，兩個(gè)旁側(cè)測(cè)點(diǎn)(-1000 m,5000 m,0 m)、(1000 m,5000 m,0 m)的響應(yīng)基本趨于背景值，而中心測(cè)點(diǎn)(0 m,5000 m,0 m)處的Ex響應(yīng)出現(xiàn)了減小，廣域視電阻率曲線也呈現(xiàn)下降.這是由兩個(gè)旁測(cè)點(diǎn)距低阻異常體較遠(yuǎn)，二次場(chǎng)響應(yīng)較弱，而中心測(cè)點(diǎn)在低阻異常體正上方導(dǎo)致的.

圖15 測(cè)點(diǎn)為(-1000 m,5000 m,0 m),(0 m,5000 m,0 m),(1000 m,5000 m,0 m)的測(cè)深曲線(a) 不同測(cè)點(diǎn)的Ex分量對(duì)比圖; (b) 不同測(cè)點(diǎn)的廣域視電阻率對(duì)比圖.Fig.15 The sounding curves at coordinates (-1000 m,5000 m,0 m), (0 m,5000 m,0 m), (1000 m,5000 m,0 m)(a) Comparison of Ex components of different measuring points; (b) Comparison of wide field apparent resistivity of different measuring points.

3 結(jié)論

本文使用單極點(diǎn)模型降階算法實(shí)現(xiàn)了三維陸地多頻可控源電磁法的快速正演.采用非結(jié)構(gòu)化網(wǎng)格進(jìn)行空間離散，結(jié)合伽遼金方法得到有限元線性方程組.選擇單個(gè)最優(yōu)化極點(diǎn)，進(jìn)一步結(jié)合直接法求解器Pardiso，利用Lanczos算法快速構(gòu)建Krylov子空間的正交基.將有限元線性方程組的系數(shù)矩陣投影到Krylov子空間，得到維度較小的降階矩陣，利用降階系統(tǒng)實(shí)現(xiàn)了多頻可控源電磁的快速正演.

本文設(shè)計(jì)了一維層狀模型和三維塊體模型，同解析解和多種數(shù)值方法的結(jié)果進(jìn)行了對(duì)比，結(jié)果表明：在滿足精度要求的條件下，只需花費(fèi)常規(guī)正演中計(jì)算一兩個(gè)頻點(diǎn)所消耗的時(shí)間，便可實(shí)現(xiàn)幾十至上百個(gè)頻率的快速正演，大大提高了正演計(jì)算的效率，驗(yàn)證了本文算法求解多頻問(wèn)題的高效性及其良好的三維適應(yīng)性.針對(duì)寬頻的問(wèn)題(log(fmax/fmin)>4)，由于單極點(diǎn)模型降階算法需要構(gòu)建更大的Krylov子空間維度以及更高的網(wǎng)格需求，我們提出頻率分段的策略，分別針對(duì)高頻和低頻部分，同時(shí)采用兩套網(wǎng)格進(jìn)行計(jì)算，這不僅可以減少網(wǎng)格不必要的浪費(fèi)，還能加快模型降階的收斂速度，以更小的Krylov子空間維度達(dá)到更高的精度.設(shè)計(jì)了三維球體模型，結(jié)合頻率分段策略進(jìn)行正演計(jì)算，數(shù)值結(jié)果表明，頻率分段策略不僅能夠?qū)崿F(xiàn)正演的快速計(jì)算，還能很好的控制高頻和低頻的精度.值得注意的是，網(wǎng)格剖分需要考慮四面體網(wǎng)格的質(zhì)量，畸變的四面體網(wǎng)格可能會(huì)導(dǎo)致病態(tài)的矩陣方程，從而影響有理Krylov子空間近似的效果.

本文開發(fā)的基于Lanczos-模型降階算法的頻率域可控源電磁正演求解器具有快速求解三維多頻可控源電磁響應(yīng)的能力，對(duì)于研究地質(zhì)模型的三維電磁響應(yīng)特征具有十分重要的意義.