基于塊狀有理Krylov方法的大地電磁三維模型降階正演

周建美 劉文韜 魯凱亮 李 貅

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

0 引言

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

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

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

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

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

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

極點是有理Krylov方法的一個重要參數，極點的選取會影響方法的計算速度和計算精度。Güttel[52]對極點的選取算法進行了總結，提出了一種最優化極點選取方法，即通過全局搜索得到多個最優化重復極點[50]。基于該方法，本文引入一種更簡單、有效的方法獲取最優化單個重復極點[42]，結合直接法求解器[53]，可實現MT三維模型的高效正演模擬。

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

1 正演算法

1.1 控制方程離散

設電磁場隨時間的變化關系為eiωt，忽略位移電流，則MT正演對應的控制方程為[14]

?×E+iωB=0

(1)

?×μ-1B-σE=S

(2)

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

Ce+iωb=0

(3)

CTMμb-Mσe=s

(4)

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

圖1 TE極化模式(左圖中藍色箭頭)和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)中的磁場，得到關于電場的控制方程

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

(8)

e(ω)=gω(A)X

(9)

1.2 塊狀有理Krylov方法

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

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

(10)

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

(11)

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

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

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

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

1.3 MT響應

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

(12)

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

(13)

(14)

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

2 數值實例

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

2.1 均勻半空間模型

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

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

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

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

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

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

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

表3 均勻半空間模型MT響應模擬不同算法計算量及耗時對比

2.2 三維DTM1模型

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

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

圖6 DTM1模型示意圖

表4 DTM1模型塊體空間坐標和電阻率

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

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

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

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

3 討論

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

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

圖9為采用BRK方法計算的51個頻率和21個頻率的MT響應結果，這里僅給出通常基于xy和yx分量計算的視電阻率和相位。如圖9所示，51個頻率的MT視電阻率和相位曲線，相比21個頻率的MT曲線，更能夠體現曲線變化的細節特征。因此，在計算時間允許的條件下，建議選取盡可能多的頻點數，可得到更精細的響應曲線特征。

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

4 結論

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

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

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