基于BTTB矩陣的快速高精度三維磁場(chǎng)正演

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

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083 2 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083 3 中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083

0 引言

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

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

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

鑒于上述存在的問(wèn)題，如何使正演算法同時(shí)具備空間域的高精度和頻率域的高效率成為新的研究熱點(diǎn).在三維重磁正反演解釋中，直立長(zhǎng)方體是最常用的一種場(chǎng)源幾何剖分單元.當(dāng)場(chǎng)源區(qū)域和觀測(cè)面都采用規(guī)則且均勻的剖分時(shí)，位場(chǎng)數(shù)據(jù)處理或正演的雅克比矩陣(即核矩陣)具有平移等效性和互換等效性(姚長(zhǎng)利等,2003)，其實(shí)質(zhì)是一類(lèi)特殊的分塊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)發(fā)展了基于BTTB矩陣的Block Circulant Extension(BCE)算法，采用FFT快速實(shí)現(xiàn)該類(lèi)矩陣和向量的卷積計(jì)算(Vogel,2002)，并用于磁場(chǎng)數(shù)據(jù)的向上和向下延拓，顯著提升了空間域數(shù)據(jù)處理的計(jì)算效率.Chen和Liu(2019)采用存儲(chǔ)中間變量的策略，進(jìn)一步優(yōu)化了BCE算法，顯著提升了重力異常正演的計(jì)算效率.Qiang等(2019)將BCE算法和分層插值法結(jié)合，分別正演了水平觀測(cè)面和起伏地形上的磁異常.Hogue等(2020)將BCE算法進(jìn)行擴(kuò)展，使其適用于磁異常正演所產(chǎn)生的非對(duì)稱(chēng)的普通BTTB矩陣.

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

1 理論和方法

1.1 基于直立長(zhǎng)方體的磁場(chǎng)正演公式

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

單個(gè)長(zhǎng)方體在任意觀測(cè)點(diǎn)P(x,y,z)產(chǎn)生的ΔT場(chǎng)及其梯度場(chǎng)的無(wú)奇點(diǎn)解析表達(dá)式為(駱遙和姚長(zhǎng)利,2007)：

(1)

(2)

(3)

(4)

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

·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場(chǎng)及其梯度場(chǎng)的核函數(shù)：

(9)

(10)

(11)

(12)

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

(13)

(14)

(15)

(16)

1.2 快速正演

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

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

(17)

(18)

(19)

(20)

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

(21)

(22)

(23)

(24)

通過(guò)預(yù)先計(jì)算并存儲(chǔ)部分中間變量可以進(jìn)一步優(yōu)化計(jì)算核矩陣的效率，引入的中間變量為：

(25)

(26)

(27)

(28)

計(jì)算一個(gè)長(zhǎng)方體單元對(duì)一個(gè)觀測(cè)點(diǎn)的磁異常時(shí)，只需要取長(zhǎng)方體單元的8個(gè)角點(diǎn)所對(duì)應(yīng)的8個(gè)中間變量進(jìn)行累加求和.當(dāng)逐層計(jì)算每一層的Nx×Ny個(gè)核矩陣元素時(shí)，如果不預(yù)先計(jì)算并存儲(chǔ)中間變量，相當(dāng)于Nz層總共計(jì)算8×Nx×Ny×Nz次中間變量.反之，則只需要計(jì)算(Nx+1)×(Ny+1)×(Nz+1)次中間變量，計(jì)算量減少至約之前的1/8.

當(dāng)Δx=Δy時(shí)，能利用對(duì)稱(chēng)性提升計(jì)算ΔT和ΔTz的核矩陣的效率.以計(jì)算ΔT為例(對(duì)ΔTz同樣成立)，一個(gè)長(zhǎng)方體單元對(duì)所有觀測(cè)點(diǎn)所產(chǎn)生的Nx×Ny個(gè)核矩陣元素具有對(duì)稱(chēng)性，如圖2所示(箭頭位置對(duì)應(yīng)的核矩陣元素相同)，用公式表示為：

圖2 一個(gè)長(zhǎng)方體單元對(duì)所有觀測(cè)點(diǎn)的磁異常ΔT核矩陣的對(duì)稱(chēng)性示意圖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)

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

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

表1對(duì)比了當(dāng)Δx=Δy，Nx=Ny時(shí)，快速正演算法分別在垂直磁化和傾斜磁化下計(jì)算ΔT場(chǎng)及其梯度場(chǎng)所需存儲(chǔ)的矩陣元素個(gè)數(shù)，據(jù)此可估算出計(jì)算量和內(nèi)存隨剖分個(gè)數(shù)的變化關(guān)系.

表1 Δx=Δy，Nx=Ny時(shí)快速正演計(jì)算方法所需存儲(chǔ)的矩陣元素個(gè)數(shù)對(duì)比Table 1 Comparison of the number of matrix elements of storage required in the proposed method，when Δx=Δy and Nx=Ny

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

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

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

2 數(shù)值模擬

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

2.1 垂直磁化模型

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

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

圖5 垂直磁化時(shí)球異常體的模型及其ΔT正演結(jié)果(a)球異常體在z=-600 m處切面的磁化率分布；(b)傳統(tǒng)解析解的計(jì)算結(jié)果；(c)快速算法的計(jì)算結(jié)果；(d)快速算法計(jì)算結(jié)果的絕對(duì)誤差.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 垂直磁化時(shí)傳統(tǒng)解析解，Gauss-FFT算法和快速算法的計(jì)算性能比較Table 2 Comparison of computational performance between the traditional analytical solution,Gauss-FFT algorithm,and the proposed method under perpendicular magnetization

類(lèi)似的，當(dāng)計(jì)算ΔTx時(shí)快速算法總耗時(shí)約2.50 s，其中核矩陣耗時(shí)約0.14 s，核矩陣與磁化率向量的乘積耗時(shí)約2.36 s，快速算法分別比傳統(tǒng)解析解方法、4點(diǎn)和8點(diǎn)Gauss-FFT算法快約3.05×104倍、26倍和105倍.圖6為傳統(tǒng)解析解和快速算法的ΔTx正演結(jié)果以及快速算法正演結(jié)果的絕對(duì)誤差，最大絕對(duì)誤差僅約5.03×10-9nT/m.

圖6 垂直磁化時(shí)球異常體的ΔTx的正演結(jié)果和誤差(a)傳統(tǒng)解析解的計(jì)算結(jié)果；(b)快速算法的計(jì)算結(jié)果；(c)快速算法計(jì)算結(jié)果的絕對(duì)誤差.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 傾斜磁化模型

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

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

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

圖7 傾斜磁化時(shí)球異常體的ΔT和ΔTx的正演結(jié)果和誤差(a)(d)、(b)(e)、(c)(f)分別為ΔT和ΔTx的傳統(tǒng)解析解的結(jié)果，快速算法的結(jié)果及其絕對(duì)誤差.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分別對(duì)比了傳統(tǒng)解析解，8點(diǎn)Gauss-FFT算法，和快速算法在計(jì)算ΔT和ΔTx時(shí)的相對(duì)計(jì)算時(shí)間隨剖分個(gè)數(shù)的變化關(guān)系，均采用最長(zhǎng)的計(jì)算時(shí)間進(jìn)行了歸一化處理.隨著剖分個(gè)數(shù)增加，正演耗時(shí)呈指數(shù)增加，而且水平方向上的剖分個(gè)數(shù)越多，快速算法相比于傳統(tǒng)解析方法提升的計(jì)算效率越高.從圖8a可知，當(dāng)計(jì)算垂直磁化下的ΔT且各方向剖分個(gè)數(shù)為240時(shí)，快速算法比傳統(tǒng)解析解方法快約5個(gè)數(shù)量級(jí)，比8點(diǎn)Gauss-FFT算法快約兩個(gè)數(shù)量級(jí).此外，垂直磁化下的快速算法比傾斜磁化下的快約一個(gè)數(shù)量級(jí).圖8b具有類(lèi)似的性質(zhì)，這也驗(yàn)證了快速算法的高效.

圖8 快速算法與傳統(tǒng)解析解、Gauss-FFT的相對(duì)計(jì)算時(shí)間對(duì)比(a)和(b)分別為ΔT和ΔTx的相對(duì)計(jì)算時(shí)間對(duì)比，都采用最長(zhǎng)的計(jì)算時(shí)間進(jìn)行了歸一化處理.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分別為垂直磁化和傾斜磁化時(shí)快速正演ΔT和ΔTx的總耗時(shí)，計(jì)算核矩陣的耗時(shí)，和核矩陣與磁化率向量的乘積耗時(shí)隨剖分個(gè)數(shù)的變化關(guān)系.從圖9a、b可知，在垂直磁化下快速正演ΔT，計(jì)算核矩陣的耗時(shí)比核矩陣與磁化率向量的乘積耗時(shí)低約一個(gè)數(shù)量級(jí)，而傾斜磁化時(shí)相反.從圖9c、d可知，快速正演ΔTx時(shí)也存在類(lèi)似現(xiàn)象.值得注意的是，垂直磁化和傾斜磁化下，快速正演中計(jì)算核矩陣與磁化率向量的乘積耗時(shí)幾乎一樣，而計(jì)算核矩陣耗時(shí)存在明顯差異.這種差異既源于垂直磁化時(shí)特殊的BTTB矩陣的結(jié)構(gòu)特性，也得益于中間變量的引入.可以看出，快速算法通過(guò)改進(jìn)核矩陣的計(jì)算過(guò)程，將正演總耗時(shí)壓縮到跟核矩陣與磁化率向量的乘積耗時(shí)相同水平.

圖9 快速算法計(jì)算ΔT和ΔTx的絕對(duì)計(jì)算時(shí)間(a)(c)、(b)(d)分別對(duì)應(yīng)ΔT和ΔTx在垂直磁化和傾斜磁化下的絕對(duì)計(jì)算時(shí)間.圖(c)的計(jì)算核矩陣耗時(shí)曲線，當(dāng)各方向剖分個(gè)數(shù)為80時(shí)實(shí)測(cè)耗時(shí)為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 應(yīng)用

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

為了驗(yàn)證快速正演算法的實(shí)用性，本文將提出的快速正演運(yùn)用于磁異常數(shù)據(jù)的聚焦反演，并對(duì)比了傳統(tǒng)聚焦反演方法與加入了快速正演算法的聚焦反演方法在計(jì)算效率和內(nèi)存需求方面的差異.

針對(duì)位場(chǎng)反演的非唯一性，Tikhonov和Arsenin(1977)提出的正則化反演目標(biāo)函數(shù)為：

φ=φd+βφm，

(30)

其中φd為數(shù)據(jù)目標(biāo)函數(shù)，定義為：

(31)

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

(32)

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

(33)

Ws為模型加權(quán)矩陣，可表示為：

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

(34)

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

利用上述公式，可得轉(zhuǎn)換后的反演目標(biāo)函數(shù)：

(35)

將公式(35)右邊第一項(xiàng)展開(kāi)得到：

=[(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)

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

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

(37)

3.2 反演模型試驗(yàn)

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

圖10 真實(shí)模型磁化率分布及其產(chǎn)生的磁異常ΔT(a)z=50 m切面磁化率分布；(b)y=300 m斷面磁化率分布；(c)加入了高斯噪聲的觀測(cè)數(shù)據(jù).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.

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

圖11 模型反演結(jié)果(a)z=50 m剖面上反演結(jié)果；(b)y=300 m斷面上反演結(jié)果.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.

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

4 結(jié)論

針對(duì)傳統(tǒng)空間域三維磁場(chǎng)正演計(jì)算效率較低的問(wèn)題，本文改進(jìn)了空間域三維磁場(chǎng)正演算法，使其具備高精度和高效率的優(yōu)點(diǎn)，并將其運(yùn)用到磁場(chǎng)反演中.基于均勻網(wǎng)格剖分，該正演方法先利用BTTB矩陣和等效性壓縮核矩陣，通過(guò)引入中間變量減少重復(fù)計(jì)算，再基于BTTB矩陣和FFT的快速算法進(jìn)一步提高計(jì)算效率.得到以下結(jié)論：

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

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

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

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

附錄A

將水平觀測(cè)面均勻剖分成nx×ny個(gè)網(wǎng)格，將場(chǎng)源區(qū)域均勻剖分成nx×ny×nz個(gè)長(zhǎng)方體單元，即垂直方向剖分成nz層，每一層在水平面又剖分成nx×ny個(gè)長(zhǎng)方體單元.假設(shè)磁化率為1，此時(shí)磁異常值只與觀測(cè)點(diǎn)與長(zhǎng)方體單元的位置有關(guān)，計(jì)算第一層單元體產(chǎn)生的磁異常值.

第一層的第一個(gè)長(zhǎng)方體單元對(duì)所有觀測(cè)點(diǎn)的磁異常用T(1,1)表示：

(A1)

其中，元素的上標(biāo)為長(zhǎng)方體單元坐標(biāo)，下標(biāo)為觀測(cè)點(diǎn)坐標(biāo)，i=1,…,nx，j=1,…,ny.僅當(dāng)nx=ny時(shí)，T(1,1)才為方陣.

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

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

(A2)

以此類(lèi)推，依次求出第一層剩下的所有長(zhǎng)方體單元對(duì)所有觀測(cè)點(diǎn)的磁異常，同樣用列向量表示.

第一層所有長(zhǎng)方體單元對(duì)所有觀測(cè)點(diǎn)的磁異常值用T1表示：

(A3)

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

(A4)

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

設(shè)真實(shí)磁化率為

(A5)

令：

(A6)

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

(A7)

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

(A8)

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

首先將t在復(fù)數(shù)域擴(kuò)展為：

(A9)

令：

(A10)

(A11)

令：

(A12)

再將磁化率進(jìn)行擴(kuò)展：

(A13)

令g=Cextvec(fext)，運(yùn)用FFT計(jì)算，其中：

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

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