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

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

0 引言

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

陸地可控源電磁的主流正演方法有積分方程法(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).有限單元法近年來得到了廣泛的應(yīng)用，主要因為它具有理論體系成熟、網(wǎng)格剖分靈活等特點，結(jié)合非結(jié)構(gòu)化網(wǎng)格可以模擬復(fù)雜地形及地質(zhì)結(jié)構(gòu)，便于在場源、測點及電性分界面處進行網(wǎng)格加密，從而提高數(shù)值計算的精度.有限單元法正演可以采用總場公式和二次場公式.總場公式需要對場源附近的電磁場做精細處理，并且在計算邊界上需要施加精細的邊界條件，但總場公式能夠有效模擬任意起伏地形(李建慧等, 2016; 邱長凱等, 2018).二次場公式的模擬變量為散射場，散射場在源附近一般不再具有奇異性，因此具有模擬精度高的特點，但是二次場公式需要背景模型具有解析的一次場表達式(張林成等, 2017).一般來說, 地形背景模型不具備解析的一次場表達式，二次場公式處理地形存在困難.因此，本文選擇總場公式作為開發(fā)可處理任意起伏地形情況的陸地可控源電磁快速正演方法.

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

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

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

1 方法

1.1 CSEM有限元線性方程組

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

(1)

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

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

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

(2)

(3)

(4)

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

(5)

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

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

(6)

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

(7)

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

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

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

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

(8)

為解決此類形如f(A)X的問題，構(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)是容易求解的，對于多頻問題，求解式(12)是快速的.

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

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)的近似表達式(12).相對Arnoldi算法而言，Lanczos算法中大型矩陣向量的乘積運算次數(shù)大大減少，可以提高構(gòu)建正交基的效率，因此本文選擇Lanczos算法來構(gòu)建正交基，算法的詳細過程在表1中給出.

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

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

(14)

2 數(shù)值實驗

本文進行數(shù)值模擬的計算環(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個人PC機.

2.1 加速比測試

2.1.1 一維模型

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

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

表2 運算時間對比Table 2 Comparison of calculation time

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

表3 3DMOR的運算時間分布Table 3 Calculation time distribution of 3DMOR

圖2 3DMOR_Lanczos同1D solution、3DFEM、3DMOR_Arnoldi的測深曲線對比(a) 不同算法的Ex響應(yīng)對比圖; (b) 3DMOR_Lanczos同其他方法之間的相對誤差.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時, 3DMOR_Lanczos同1D solution、3DFEM、3DMOR_Arnoldi對比(a) 不同算法的Ex響應(yīng)對比圖; (b) 3DMOR_Lanczos同其他方法之間的相對誤差.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 三維模型

為進一步說明本文算法的高效性和對三維模型的適應(yīng)性，設(shè)計如圖4所示的三維模型進行對比測試.圖4a為模型的剖面圖，圖4b為模型的平面圖，圖中虛線為測線.在電阻率為50 Ωm的均勻半空間中嵌入一個大小為120 m×200 m×400 m的低阻異常體，其中心坐標為(1000 m,0 m,300 m)，電阻率為5 Ωm.100 m長的源沿著x軸放置，源的中心坐標為(50 m,0 m,0 m)，發(fā)射電流為1 A，發(fā)射頻率為1～10000 Hz對數(shù)等間隔分布的32個頻點.網(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)、基于磁矢勢的有限元方法(FE)(Ansari and Farquharson, 2014)、積分方程法(IE)(Farquharson and Oldenburg, 2002; Zhou et al., 2018a)的數(shù)值解進行了Ex實部和虛部響應(yīng)曲線的對比，正演計算頻率為3 Hz，如圖5a、圖6a所示.從圖中可知，由于地下低阻體的存在，地下電流被低阻體“吸引”，導(dǎo)致地表的電流密度降低，從而出現(xiàn)了電場Ex響應(yīng)的實部和虛部曲線在x為900～1100 m范圍內(nèi)均呈現(xiàn)向下凹陷現(xiàn)象.圖5b和圖6b給出了Ex響應(yīng)的相對誤差曲線，除IE方法外，本文算法同其他兩種方法之間的相對誤差未超過3%，這說明利用本文算法進行三維模型的正演是可行的.

圖5 頻率為3 Hz時, 3DMOR同多種數(shù)值模擬方法的電場實部的對比(a) 不同算法的Ex實分量對比圖; (b) 3DMOR同其他方法之間的相對誤差.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時, 3DMOR同多種數(shù)值模擬方法的電場虛部的對比(a) 不同算法的Ex虛分量對比圖; (b) 3DMOR同其他方法之間的相對誤差.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給出了測點(1000 m,0 m,0 m)處3DMOR和3DFEM的Ex響應(yīng)曲線及相對誤差曲線，可以看出3DMOR對3DFEM的近似效果較好.同樣在表4給出了計算資源消耗的對比，從中可知加速比達到了17.6.與表2相比可知，隨著未知數(shù)的增加，3DMOR更節(jié)省內(nèi)存.通過表3我們可以發(fā)現(xiàn)，隨著未知數(shù)的增加，矩陣分解和正交化過程將消耗大部分的運算時間，因此本文采用的高效的Lanczos正交化算法是可以提升正演計算效率的.

圖7 中心測點的測深曲線，3DMOR同3DFEM的對比(a) 兩種算法的Ex響應(yīng)對比圖; (b) 3DMOR同3DFEM之間的相對誤差.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 運算時間對比Table 4 Comparison of calculation time

2.2 頻率分段策略

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

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

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

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

圖9 當頻率計算范圍為[0.01,10000]Hz時，3DMOR采取不同Krylov子空間維度的測深曲線同1D solution對比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

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

圖10 當頻率計算范圍為[0.01,10000]Hz時，采取頻率分段策略的計算結(jié)果(a) 3DMOR同1D solution的測深曲線對比; (b) 相對誤差.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è)計如圖11所示的三維低阻球體模型，圖12為網(wǎng)格剖分示意圖.球體的半徑為200 m，中心埋深位于(0 m,5000 m,400 m)，電阻率為10 Ωm，圍巖電阻率為100 Ωm.x方向的電流源放置在坐標原點，其長度為200 m，發(fā)射電流為10 A，發(fā)射頻率為0.01～10000 Hz對數(shù)等間隔分布的100個頻點.本例在地表布置了11條測線，線距為200 m，y坐標范圍為-1000～1000 m.每條測線長度為2000 m，含11個測點，點距為200 m，x坐標范圍為4000～6000 m.仍然采用上文的頻率分段策略來進行正演求解.對于高頻段，設(shè)置網(wǎng)格邊界為±20 km，方程未知數(shù)為578603，構(gòu)建Krylov子空間的維度為80.對于低頻段，網(wǎng)格邊界為±150 km，方程未知數(shù)為387979，構(gòu)建Krylov子空間的維度為80.兩套網(wǎng)格并行計算，最終消耗總內(nèi)存為4.11 GB，總運行時間為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給出了四個頻率的磁場Hy分量、電場Ex分量以及廣域視電阻率的切片圖，四個頻率分別為10000 Hz、100 Hz、1 Hz、0.01 Hz，其中廣域視電阻率利用Ex求得(李帝銓等, 2013; He, 2018;武建平等, 2020；Yu et al., 2020).從圖13(a,b,c,d)中可以看出，磁場Hy分量對異常體的響應(yīng)較微弱，無法從切片圖中觀察到異常場造成的畸變.如圖13(f,g,h,j,k,l)所示，由于低阻異常體的存在，使得地表電流密度減小，因此電場Ex分量和廣域視電阻率均出現(xiàn)凹陷，在頻率為100 Hz、1 Hz、0.01 Hz均有不同程度的響應(yīng)，能夠較好的反映異常體的位置和屬性.圖14給出的是兩條中心測線(x=0 m測線,y=5000 m測線)的廣域視電阻率擬斷面圖，可以看出，在高頻部分廣域視電阻率趨向于背景電阻率，在低頻部分呈現(xiàn)低阻響應(yīng)形態(tài).

圖13 頻率分別為10000 Hz、100 Hz、1 Hz、0.01 Hz的Hy、Ex、廣域視電阻率切片圖(a)—(d)磁場Hy分量切片圖; (e)—(h)電場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測線和y=5000測線的廣域視電阻率擬斷面圖Fig.14 Pseudo-section of wide field apparent resistivity with line x=0 and line y=5000

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

圖15 測點為(-1000 m,5000 m,0 m),(0 m,5000 m,0 m),(1000 m,5000 m,0 m)的測深曲線(a) 不同測點的Ex分量對比圖; (b) 不同測點的廣域視電阻率對比圖.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é)論

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

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

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