基于BTTB矩陣的快速高精度三維磁場正演

袁洋，崔益安，陳波*，趙廣東，柳建新，郭榮文

1 中南大學地球科學與信息物理學院，長沙 410083 2 有色資源與地質災害探查湖南省重點實驗室，長沙 410083 3 中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室，長沙 410083

0 引言

磁法勘探是一種最為常見的地球物理勘探方法，其主要利用地下巖石和礦物的磁性差異所引起的磁異常來區分目標異常體(柳建新等,2016)，已被廣泛應用于礦產資源勘探(Hinze et al.,2013；周文月等，2021)、考古(Cella and Fedi,2015;Fedi and Florio,2003)、軍事地球物理(Hiergeist et al.,2015)、區域大地構造和板塊碰撞演化(Gao et al.,2013;Hemant and Mitchell,2009;Rajaram and Langel,1992；羅凡等，2021)等方面.隨著高精度衛星磁測和航空磁測的普及，磁法勘探將會在越來越多的領域中發揮重要作用.通過三維磁化率成像反演，高精度和高分辨率的磁測數據能反映地下三維磁化率分布，并獲得豐富的地質構造信息.但當處理大規模磁異常數據時，磁異常三維反演算法存在計算效率較低和內存占用巨大等問題，嚴重制約著其在磁法勘探的實際應用.

正演是反演的基礎，正演算法的計算效率和計算精度直接影響著反演的耗時和可信度.磁場正演計算方法可分為空間域方法(Vacquier et al.,1951;Hughson,1964;Bhattacharyya,1964;Sharma,1966;Hjelt,1972;Plouff,1976;Kunaratnam,1981;Blakely,1996;Ren et al.,2017a,b;郭志宏等,2004;駱遙和姚長利,2007)和頻率域方法(Bhattacharyya,1966;Tontini,2005;Tontini et al.,2009;吳宣志,1981;柴玉璞和萬海珍,2020).空間域算法通常采用解析公式單獨計算各個單元體的異常值，再進行累加求和.這類算法具有計算精度高的優點，但其計算效率較低，特別是當模型復雜需要精細剖分時，計算時間會顯著增加.

頻率域方法由于引入了快速傅里葉變換(FFT)，計算效率較高，但受限于FFT算法固有的缺陷，例如混疊效應、邊界效應和強制周期化等(Bracewell,1978;Percival and Walden,1993;Boyd,2001;Wu and Tian,2014)，導致這類算法正演精度較低.近年來，頻率域算法也不斷得到改進和運用.例如，Wu和Tian(2014)提出高斯-快速傅立葉變換(Gauss-FFT)算法，在每一個積分區間采用高精度高斯型數值積分代替傳統FFT算法中的梯形法則，大大提高了正演計算精度，被廣泛用于頻率域三維重磁場正演(Zhao et al.,2018;曾明等,2019)、莫霍面反演(Zhao et al.,2020)和地形校正(Wu,2016;Wu and Chen,2016)等.類似地，基于FFT的空間波數域混合方法(李昆等,2019;周印明等,2020)同樣顯著提升了正演的計算效率.雖然這些頻率域方法在較高計算效率的情況下實現了較高的計算精度，但是與空間域解析解之間的誤差還是不可忽略.

鑒于上述存在的問題，如何使正演算法同時具備空間域的高精度和頻率域的高效率成為新的研究熱點.在三維重磁正反演解釋中，直立長方體是最常用的一種場源幾何剖分單元.當場源區域和觀測面都采用規則且均勻的剖分時，位場數據處理或正演的雅克比矩陣(即核矩陣)具有平移等效性和互換等效性(姚長利等,2003)，其實質是一類特殊的分塊Toeplitz矩陣，即塊-托普利茲-托普利茲-塊(Block-Toeplitz Toeplitz-Block,BTTB)矩陣(陳龍偉,2013;Zhang and Wong,2015;Zhang et al.,2016;Wu,2018;Chen and Liu,2019;Hogue et al.,2020).陳龍偉(2013)發展了基于BTTB矩陣的Block Circulant Extension(BCE)算法，采用FFT快速實現該類矩陣和向量的卷積計算(Vogel,2002)，并用于磁場數據的向上和向下延拓，顯著提升了空間域數據處理的計算效率.Chen和Liu(2019)采用存儲中間變量的策略，進一步優化了BCE算法，顯著提升了重力異常正演的計算效率.Qiang等(2019)將BCE算法和分層插值法結合，分別正演了水平觀測面和起伏地形上的磁異常.Hogue等(2020)將BCE算法進行擴展，使其適用于磁異常正演所產生的非對稱的普通BTTB矩陣.

為了進一步降低磁場正反演中核矩陣的內存占用并提高計算效率，本文借鑒Chen和Liu(2019)重力場正演的優化方法，改進了空間域磁異常ΔT及其梯度的正演算法，使其更為高效.首先采用無解析奇點的長方體ΔT場及其梯度場公式保證計算精度，利用BTTB矩陣的等效性壓縮核矩陣大大降低內存需求，并且在計算核矩陣時引入中間變量減少重復計算；然后采用Vogel(2002)中基于BTTB矩陣和FFT的快速算法進一步提高計算效率；再分別進行垂直磁化和傾斜磁化兩種情況下的正演數值模擬，驗證算法的效率和精度；最后將快速算法運用到磁場反演中，并用合成模型實驗驗證算法的實用性.

1 理論和方法

1.1 基于直立長方體的磁場正演公式

本研究的場源區域和觀測面的剖分如圖1所示，坐標原點位于場源區域的上界面，z軸向下.將地下場源區域沿x、y、z方向剖分成Nx×Ny×Nz個長方體單元，三個方向間隔分別為Δx、Δy、Δz.(ξa,ηb,ζc)為長方體單元(a,b,c)的中心坐標.假設每個長方體單元的磁化率為常數，κ(ξa,ηb,ζc)為該單元體對應的磁化率，其中a=1,2,…,Nx，b=1,2,…,Ny，c=1,2,…,Nz.觀測點分布于高度為z0平面上的Nx×Ny個規則網格節點上.

單個長方體在任意觀測點P(x,y,z)產生的ΔT場及其梯度場的無奇點解析表達式為(駱遙和姚長利,2007)：

(1)

(2)

(3)

(4)

值得注意的是，公式(1)計算的ΔT為磁異?？倧姸仁噶縏a在T0方向上的投影.對于一般情況(即當磁異常強度不大時)，實測的ΔT(即磁場總強度T與T0的模量差)可近似為式(1)正演計算的ΔT，通常具有足夠精度.但當存在強磁性體或磁異常幅值大時，實測ΔT和正演ΔT之間的誤差不可忽略(袁曉雨等，2015).一些學者提出了實測ΔT的處理和正反演的方法(Zhen et al.,2019;Sun et al.,2019;甄慧翔等,2019;孫石達等,2020;胡正旺等,2021)，有效改善了兩者之間的誤差.本文針對一般情況，根據線性疊加原理，地下場源在觀測點(xi,yj,z0)產生的ΔT場及其梯度場可以表示為所有長方體單元磁場響應累加求和的形式：

·G(xi,yj,z0;ξa,ηb,ζc)，

(5)

·Gx(xi,yj,z0;ξa,ηb,ζc)，

(6)

·Gy(xi,yj,z0;ξa,ηb,ζc)，

(7)

·Gz(xi,yj,z0;ξa,ηb,ζc)，

(8)

式中，i=1,2,…,Nx，j=1,2,…,Ny，G、Gx、Gy、Gz分別為ΔT場及其梯度場的核函數：

(9)

(10)

(11)

(12)

公式(5)—(8)可以采用矩陣形式表示為(沿場源區域的垂直方向先逐層計算再進行累加求和)：

(13)

(14)

(15)

(16)

1.2 快速正演

基于上述剖分，當觀測點與長方體單元中心點在水平面上的投影重合時，對于任意的Nx、Ny、Δx、Δy，上述核矩陣都是一類特殊的矩陣，稱為BTTB矩陣，即分塊Toeplitz矩陣.這類矩陣具有大量重復的矩陣元素，因此僅需計算部分元素即可得到全部矩陣，從而可以減少正演過程中計算核矩陣的耗時.例如，當場源為傾斜磁化時，ΔT及其梯度的核矩陣為普通的BTTB矩陣，在正演時僅需計算(2Nx-1)×(2Ny-1)個核矩陣元素就能得到完整的核矩陣.

特別是垂直磁化時，組成Gc、Gzc、Gyc的對稱Toeplitz矩陣是對稱分布的，但Gyc中對稱的兩個Toeplitz矩陣的元素互為相反數.另外，組成Gxc的Toeplitz矩陣塊也是對稱分布的，矩陣塊中沿主對角線對稱的元素互為相反數.對于這些特殊的BTTB矩陣，只需要計算核矩陣的第一行或者第一列(即Nx×Ny個核矩陣元素)就能構建完整的核矩陣，當計算第一層長方體單元對所有觀測點產生的磁異常時，只需第一層第一個單元體(ξ1,η1,ζ1)對所有觀測點(xi,yj,z0)所產生的Nx×Ny個核矩陣元素，i=1,2,…,Nx，j=1,2,…,Ny，即令：

(17)

(18)

(19)

(20)

對于垂直磁化(I0=I=90°、A0=A=0°)，此時l=m=L=M=0、n=N=1，于是K1=K2=K4=K5=K6=0、K3=2，公式(17)—(20)可簡化成：

(21)

(22)

(23)

(24)

通過預先計算并存儲部分中間變量可以進一步優化計算核矩陣的效率，引入的中間變量為：

(25)

(26)

(27)

(28)

計算一個長方體單元對一個觀測點的磁異常時，只需要取長方體單元的8個角點所對應的8個中間變量進行累加求和.當逐層計算每一層的Nx×Ny個核矩陣元素時，如果不預先計算并存儲中間變量，相當于Nz層總共計算8×Nx×Ny×Nz次中間變量.反之，則只需要計算(Nx+1)×(Ny+1)×(Nz+1)次中間變量，計算量減少至約之前的1/8.

當Δx=Δy時，能利用對稱性提升計算ΔT和ΔTz的核矩陣的效率.以計算ΔT為例(對ΔTz同樣成立)，一個長方體單元對所有觀測點所產生的Nx×Ny個核矩陣元素具有對稱性，如圖2所示(箭頭位置對應的核矩陣元素相同)，用公式表示為：

圖2 一個長方體單元對所有觀測點的磁異常ΔT核矩陣的對稱性示意圖Fig.2 The schematic diagram of the symmetry of the kernel matrix generated by a cuboid unit for the magnetic anomaly ΔT

G(xi,yj,z0;ξ1,η1,ζ1)=G(xj,yi,z0;ξ1,η1,ζ1).

(29)

事實上，正演計算并不需要一次性計算并存儲整個Nz+1個界面的中間變量，而是在逐層計算每一層核矩陣元素時，僅預先計算并存儲這一層的中間變量，對應長方體單元上下兩個界面的角點.這樣增加的內存與三維磁化率向量相比數量級較小，可忽略不計.圖3直觀地展示了中間變量的功能，每個長方體單元在其上、下界面各有4個角點，而且上下兩個長方體單元會共用4個角點，例如長方體單元(ξ1,η1,ζ1)與(ξ1,η1,ζ2)共用4個角點(白色圓點表示).為了說明中間變量的互換性，以長方體單元(ξ1,η1,ζ1)的上界面4個角點(三角形表示)為例，當觀測點位于對角線上的網格時，關于對角線對稱的角點對應的t值相同，即實線對應的t值相同(實線，虛線和點劃線分別對應3個不同的t值).

圖3 中間變量t值等效示意圖Fig.3 The schematic diagram of the equivalence of t values

表1對比了當Δx=Δy，Nx=Ny時，快速正演算法分別在垂直磁化和傾斜磁化下計算ΔT場及其梯度場所需存儲的矩陣元素個數，據此可估算出計算量和內存隨剖分個數的變化關系.

表1 Δx=Δy，Nx=Ny時快速正演計算方法所需存儲的矩陣元素個數對比Table 1 Comparison of the number of matrix elements of storage required in the proposed method，when Δx=Δy and Nx=Ny

利用上述過程構建磁異常ΔT及其梯度的正演核矩陣后，可以利用基于BTTB矩陣和FFT的快速算法(Vogel，2002)計算該層核矩陣與其對應的磁化率向量的乘積，進一步提高正演計算效率，BTTB矩陣和向量乘積的快速算法的詳細過程見附錄A.

圖4為快速正演算法流程圖，其中虛線路線提升的計算效率最高，點劃線最低.

圖4 快速正演算法流程圖Fig.4 The flow chart of the fast forward method

2 數值模擬

為了驗證快速算法的準確性和高效性，本節將從計算效率，計算精度和內存需求這三個方面對傳統空間域解析解(Blakely,1996)，最新頻率域Gauss-FFT算法(Wu and Tian,2014)和改進的快速算法進行對比.由于使用的高斯點個數不同，Gauss-FFT算法的計算效率和計算精度也不同，為了更合理的對比計算效率，本文將選用4點和8點的Gauss-FFT算法.另外，從1.2節可知，相比于傾斜磁化，垂直磁化時形成的BTTB矩陣具有更多重復的元素，所用的快速算法計算效率更高，因此本節會分開進行討論.

2.1 垂直磁化模型

本文以計算球異常體在水平觀測面產生的磁異常ΔT及其x方向梯度ΔTx為例進行說明.采用垂直磁化，場源區域(0 m,0 m,0 m)～(1200 m,1200 m,1200 m)，剖分個數為240 × 240 × 240，網格間距為5 m × 5 m，觀測點個數為240 × 240，觀測面高度z0=-10 m，z軸垂直向下為正.如圖5a所示，球異常體的中心位置(600 m,600 m,600 m)，半徑為200 m，異常體磁化率κ=0.03 SI，背景磁化率為0 SI，正常場大小T0=50000 nT.計算機參數為：AMD Ryzen 7 3700X 3.59 GHz(CPU)，內存32 GB.

表2為垂直磁化時傳統解析解，Gauss-FFT算法和快速算法的計算性能比較.快速算法計算ΔT總耗時約2.67 s，其中核矩陣耗時約0.29 s，核矩陣與磁化率向量的乘積耗時約2.38 s，內存需求約111.29 MB(其中，三維磁化率向量約110.59 MB，核矩陣約0.23 MB，中間變量約0.47 MB)，快速算法的內存組成部分可參考表1.相同的模型，快速算法的計算效率比傳統解析解方法快約1.55×105倍，比4點Gauss-FFT算法快約23倍，8點Gauss-FFT算法時擴大到了94倍，而且快速算法的計算精度比4點和8點的Gauss-FFT算法更高，其中4點Gauss-FFT算法的最大絕對誤差比其他三種算法大幾個數量級.以上可知，本文改進的快速算法在計算性能方面具有顯著優勢.此外，由表1可知，在開展磁異常反演時，若采用傳統解析解方法進行正演，則反演時需存儲的核矩陣約6.37×106MB，而采用快速算法僅約55.53 MB，內存需求降低了約1.15×105倍.圖5b、c分別為傳統解析解和快速算法的ΔT正演結果，圖5d為快速算法正演結果的絕對誤差，最大絕對誤差僅約1.01×10-6nT.

圖5 垂直磁化時球異常體的模型及其ΔT正演結果(a)球異常體在z=-600 m處切面的磁化率分布；(b)傳統解析解的計算結果；(c)快速算法的計算結果；(d)快速算法計算結果的絕對誤差.Fig.5 A spherical anomalous model and the ΔT results and the errors under perpendicular magnetization(a)The susceptibility distribution in the section at z=-600 m；(b)—(d)are the ΔT results of the traditional analytical solution,the proposed method,and the absolute errors.

表2 垂直磁化時傳統解析解，Gauss-FFT算法和快速算法的計算性能比較Table 2 Comparison of computational performance between the traditional analytical solution,Gauss-FFT algorithm,and the proposed method under perpendicular magnetization

類似的，當計算ΔTx時快速算法總耗時約2.50 s，其中核矩陣耗時約0.14 s，核矩陣與磁化率向量的乘積耗時約2.36 s，快速算法分別比傳統解析解方法、4點和8點Gauss-FFT算法快約3.05×104倍、26倍和105倍.圖6為傳統解析解和快速算法的ΔTx正演結果以及快速算法正演結果的絕對誤差，最大絕對誤差僅約5.03×10-9nT/m.

圖6 垂直磁化時球異常體的ΔTx的正演結果和誤差(a)傳統解析解的計算結果；(b)快速算法的計算結果；(c)快速算法計算結果的絕對誤差.Fig.6 The forward modeling results (ΔTx)and the error for the spherical anomalous body under perpendicular magnetization(a)—(c)are the ΔTx results of the traditional analytical solution,the proposed method and the absolute error.

2.2 傾斜磁化模型

為了驗證快速算法在傾斜磁化情況下的計算性能，仍采用上述球異常體模型，計算該模型在傾斜磁化下(I0=I=45°、A0=A=5°)的ΔT和ΔTx.

表3為傾斜磁化時傳統解析解，Gauss-FFT算法和快速算法的計算性能比較.快速算法計算ΔT總耗時約31.89 s，其中核矩陣耗時約29.19 s，計算ΔTx總耗時約8.45 s，其中核矩陣耗時約5.77 s.對比表2，在垂直磁化和傾斜磁化兩種情況下，傳統解析解方法和8點Gauss-FFT算法的計算性能幾乎不變.計算ΔT時快速算法分別比傳統解析解方法和8點Gauss-FFT算法快約1.33×104倍和7倍，計算ΔTx時則快約1.02×104倍和31倍.由于高斯點數的減少，4點Gauss-FFT算法比8點Gauss-FFT算法快約3倍，但計算精度較低.圖7分別為傾斜磁化時ΔT和ΔTx的正演計算結果和誤差，可知快速算法的計算精度很高.

表3 傾斜磁化時傳統解析解，Gauss-FFT算法，和快速算法的計算性能比較Table 3 Comparison of computational performance between the traditional analytical solution,Gauss-FFT algorithm and the proposed method under oblique magnetization

圖7 傾斜磁化時球異常體的ΔT和ΔTx的正演結果和誤差(a)(d)、(b)(e)、(c)(f)分別為ΔT和ΔTx的傳統解析解的結果，快速算法的結果及其絕對誤差.Fig.7 The forward modeling results (ΔT and ΔTx)and their errors for a spherical anomalous body under oblique magnetization(a)(d),(b)(e),(c)(f)are the ΔT and ΔTx results of the traditional analytical solution,the proposed method and the absolute errors.

圖8分別對比了傳統解析解，8點Gauss-FFT算法，和快速算法在計算ΔT和ΔTx時的相對計算時間隨剖分個數的變化關系，均采用最長的計算時間進行了歸一化處理.隨著剖分個數增加，正演耗時呈指數增加，而且水平方向上的剖分個數越多，快速算法相比于傳統解析方法提升的計算效率越高.從圖8a可知，當計算垂直磁化下的ΔT且各方向剖分個數為240時，快速算法比傳統解析解方法快約5個數量級，比8點Gauss-FFT算法快約兩個數量級.此外，垂直磁化下的快速算法比傾斜磁化下的快約一個數量級.圖8b具有類似的性質，這也驗證了快速算法的高效.

圖8 快速算法與傳統解析解、Gauss-FFT的相對計算時間對比(a)和(b)分別為ΔT和ΔTx的相對計算時間對比，都采用最長的計算時間進行了歸一化處理.Fig.8 Comparison of the relative computation time between the traditional analytical solution,Gauss-FFT algorithm,and the proposed methodThe time of ΔT (a)and ΔTx (b)is normalized by the longest computation time.

圖9分別為垂直磁化和傾斜磁化時快速正演ΔT和ΔTx的總耗時，計算核矩陣的耗時，和核矩陣與磁化率向量的乘積耗時隨剖分個數的變化關系.從圖9a、b可知，在垂直磁化下快速正演ΔT，計算核矩陣的耗時比核矩陣與磁化率向量的乘積耗時低約一個數量級，而傾斜磁化時相反.從圖9c、d可知，快速正演ΔTx時也存在類似現象.值得注意的是，垂直磁化和傾斜磁化下，快速正演中計算核矩陣與磁化率向量的乘積耗時幾乎一樣，而計算核矩陣耗時存在明顯差異.這種差異既源于垂直磁化時特殊的BTTB矩陣的結構特性，也得益于中間變量的引入.可以看出，快速算法通過改進核矩陣的計算過程，將正演總耗時壓縮到跟核矩陣與磁化率向量的乘積耗時相同水平.

圖9 快速算法計算ΔT和ΔTx的絕對計算時間(a)(c)、(b)(d)分別對應ΔT和ΔTx在垂直磁化和傾斜磁化下的絕對計算時間.圖(c)的計算核矩陣耗時曲線，當各方向剖分個數為80時實測耗時為0，故缺省.Fig.9 The absolute computation time of ΔT and ΔTx by the proposed method(a)(c),(b)(d)represent ΔT and ΔTx under perpendicular magnetization and oblique magnetization,respectively.As for the time-consuming curve of calculating nuclear matrix in figure (c),when the number of divisions in each direction is 80,the measured time-consuming is 0,so it is defaulted.

3 應用

3.1 基于快速正演的聚焦反演

為了驗證快速正演算法的實用性，本文將提出的快速正演運用于磁異常數據的聚焦反演，并對比了傳統聚焦反演方法與加入了快速正演算法的聚焦反演方法在計算效率和內存需求方面的差異.

針對位場反演的非唯一性，Tikhonov和Arsenin(1977)提出的正則化反演目標函數為：

φ=φd+βφm，

(30)

其中φd為數據目標函數，定義為：

(31)

其中Wd為數據加權矩陣，Wd=diag[1/ε]，其中εi為第i個觀測數據誤差的標準差.dobs為磁異常觀測數據，β為正則化參數.φm為模型目標函數，Portniaguine和Zhdanov(1999)提出了最小支撐穩定函數用來進行聚焦反演.為了減少聚焦反演迭代次數，本文采用Rezaie(2020)改進后的sigmoid穩定函數進行反演：

(32)

其中σ為聚焦參數，一般取接近于0的小數.模型目標函數可以用L2范數表示為：

(33)

Ws為模型加權矩陣，可表示為：

Ws=diag[(1+exp(-κ2/σ2))(κ2+σ2)]-0.5，

(34)

Wz為深度加權矩陣，Wz=diag[1/z1.5]，其中z為模型參數的深度(Li and Oldenburg,2003).

利用上述公式，可得轉換后的反演目標函數：

(35)

將公式(35)右邊第一項展開得到：

=[(Gwκw)T-(dw)T](Gwκw-dw)

=[(κw)T(Gw)T-(dw)T](Gwκw-dw)

=(κw)T[(Gw)T(Gwκw-dw)]

-(dw)T(Gwκw-dw)，

(36)

其中最消耗內存和計算時間的部分，即Gwκw中核矩陣和磁化率向量的乘積運算Gκ，以及(Gw)T中核矩陣的轉置GT和Wd的乘積運算，都引入快速正演方法優化了核矩陣內存需求，并提升了計算效率.

為了使目標函數公式(35)最小，Portniaguine和Zhdanov(1999)提出了重加權正則化共軛梯度法(Reweighted Regularized Conjugate Gradient，RRCG)計算最優解κw.本文采用RRCG算法的進行反演，詳細實現過程可參考Rezaie(2020)，同樣加入了自適應正則化方法并強制約束模型參數的上下邊界.停止迭代的判斷條件為：

(37)

3.2 反演模型試驗

本節采用合成模型進行測試，反演垂直磁化下的磁異常ΔT.場源區域(0 m,0 m,0 m)～(600 m,600 m,300 m)，剖分個數120×120×60，網格間距5 m×5 m，觀測點個數120×120，觀測面高度z0=-0.1 m，z軸垂直向下為正.階梯異常體模型在z=50 m和y=300 m處的磁化率分布分別如圖10a、b所示.首先，利用上述基于BTTB矩陣的快速正演算法，計算該模型在觀測面上產生的理論磁異常ΔT，并加入均值為零，標準差為5%理論值最大值的高斯隨機噪聲，作為反演的觀測數據，磁異常如圖10c所示.

圖10 真實模型磁化率分布及其產生的磁異常ΔT(a)z=50 m切面磁化率分布；(b)y=300 m斷面磁化率分布；(c)加入了高斯噪聲的觀測數據.Fig.10 The synthetic model and its ΔT result(a)Magnetic susceptibility distribution on horizontal section z=50 m;(b)Magnetic susceptibility distribution on longitudinal cross profile x=300 m;(c)Observation data adding Gaussian random noises.

反演計算中，采用的聚焦參數σ為0.005，磁化率約束范圍為0 ～ 0.06 SI，自適應正則化的衰減系數為0.6，收斂閾值δ為10-3.最終反演迭代了1018次，反演結果在z=50 m和y=300 m處的磁化率分布分別如圖11a、b所示.可以看出，使用RRCG算法的聚焦反演能較好的擬合真實模型的磁化率，頂面，底面和走向.另外，測試發現，聚焦參數和強制約束模型參數的上邊界的選取都對反演結果影響較大.

圖11 模型反演結果(a)z=50 m剖面上反演結果；(b)y=300 m斷面上反演結果.Fig.11 Inversed magnetic susceptibility distribution of the model(a)Inversed result in z=50 m section;(b)Inversed result in y=300 m profile.

為了對比傳統反演方法與快速反演方法的計算效率，分別測試兩者在反演過程中迭代一次的耗時，傳統反演耗時約568.13 s，快速反演耗時約0.84 s，減少了6.75×102倍.另外，傳統反演的核矩陣內存需求約為9.95×104MB，快速反演則約為3.48 MB，降低了約2.88×104倍.為了能進一步對比更大的模型，可利用等效性壓縮核矩陣的內存，解除物理內存對傳統反演方法的使用限制.當剖分個數為240×240×120，觀測點個數為240×240時，在其他反演參數不變的情況下，反演結果幾乎不變.此時，迭代一次，傳統反演耗時約2.36×104s，快速反演耗時約7.42 s，相差3.18×103倍.

4 結論

針對傳統空間域三維磁場正演計算效率較低的問題，本文改進了空間域三維磁場正演算法，使其具備高精度和高效率的優點，并將其運用到磁場反演中.基于均勻網格剖分，該正演方法先利用BTTB矩陣和等效性壓縮核矩陣，通過引入中間變量減少重復計算，再基于BTTB矩陣和FFT的快速算法進一步提高計算效率.得到以下結論：

(1)正演球異常體模型的實驗表明，計算ΔT時，垂直磁化下快速算法比傳統解析解快約1.55×105倍，比8點Gauss-FFT算法快約94倍，最大絕對誤差約為1.01×10-6nT；傾斜磁化下則分別快約1.33×104倍和7倍，最大絕對誤差約為5.54×10-9nT.證明了快速算法具有計算效率高、精度高的優點.另外，快速算法在垂直磁化時提升的計算效率比傾斜磁化時更高，說明改進核矩陣的計算過程對于提升計算效率有很大幫助.

(2)將快速正演算法引入到反演方法中，并且利用快速正演優化了反演中有核矩陣參與計算的部分.合成模型的反演結果表明，加入了快速正演算法的反演方法比傳統反演方法快了約6.75×102倍.同時，核矩陣的內存需求降低了約2.88×104倍.這也證明了快速算法具有內存需求小的優點.

(3)事實上，計算時間和內存需求最終取決于剖分個數.隨著剖分個數增加，正演耗時呈指數增加，而且水平方向上的剖分個數越多，快速算法比傳統解析解方法提升的計算效率越高.相應的，反演時壓縮的核矩陣內存越多.本文算例均采用串行計算，由于在深度方向上為逐層計算再累加求和，計算高度并行，若采用并行，計算效率還將進一步提高.鑒于快速算法在計算效率上的顯著優勢，還可將其用于大規模重力場正反演和位場延拓等處理計算，具有廣闊的應用前景.

致謝感謝三位匿名審稿人提供的建設性的意見和建議，使本文更加完善；感謝陳龍偉教授、杜勁松教授和孫石達博士在論文研究開展過程中提供的幫助和寶貴建議.

附錄A

將水平觀測面均勻剖分成nx×ny個網格，將場源區域均勻剖分成nx×ny×nz個長方體單元，即垂直方向剖分成nz層，每一層在水平面又剖分成nx×ny個長方體單元.假設磁化率為1，此時磁異常值只與觀測點與長方體單元的位置有關，計算第一層單元體產生的磁異常值.

第一層的第一個長方體單元對所有觀測點的磁異常用T(1,1)表示：

(A1)

其中，元素的上標為長方體單元坐標，下標為觀測點坐標，i=1,…,nx，j=1,…,ny.僅當nx=ny時，T(1,1)才為方陣.

將T(1,1)轉換為列向量：

T(1,1)=vec(T(1,1))

(A2)

以此類推，依次求出第一層剩下的所有長方體單元對所有觀測點的磁異常，同樣用列向量表示.

第一層所有長方體單元對所有觀測點的磁異常值用T1表示：

(A3)

其中T1為BTTB矩陣，可以表示成T1=BTTB(t)，其中：

(A4)

得到大小為(2nx-1)×(2ny-1)的矩陣t是算法非常關鍵的一步.

設真實磁化率為

(A5)

令：

(A6)

公式(A5)可用二維離散卷積表示為：

(A7)

其中，i=1,…,nx，j=1,…,ny，即：

(A8)

而二維離散卷積可以利用FFT加速算法.

首先將t在復數域擴展為：

(A9)

令：

(A10)

(A11)

令：

(A12)

再將磁化率進行擴展：

(A13)

令g=Cextvec(fext)，運用FFT計算，其中：

g=Cextvec(fext)=vec(F-1{F{cext}*F{fext}})=vec(ifft(fft2(cext).*fft2(fext))),(A14)

令g′=array(g)，抽取g′的前nx×ny個元素并排列成列向量就得到h1.類似的，將每一層長方體單元的磁異常進行疊加得到總的磁異常.