周印明 戴世坤 李 昆 凌嘉宣 胡曉穎 熊 彬

（①中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長沙410083；②中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室，湖南長沙410083；③東方地球物理公司綜合物化探處，河北涿州072751；④桂林理工大學(xué)地球科學(xué)學(xué)院，廣西桂林541006）

0 引言

重磁位場正演方法主要分為空間域方法和頻率域方法?？臻g域方法又可以分為解析法和數(shù)值法。解析法是通過重磁異常場解析表達(dá)式直接求取空間任意點(diǎn)的精確異常場，其優(yōu)點(diǎn)是原理簡單、計(jì)算結(jié)果精確。解析法研究早期側(cè)重于簡單、規(guī)則模型的重磁場解析表達(dá)式的推導(dǎo)，如水平薄板、直立棱柱體、直立圓柱體等［1-4］。后來學(xué)者們逐步對均勻介質(zhì)的多面體模型［5-9］和簡單物性變化的任意三度體模型的重磁場解析式推導(dǎo)［10-11］、解析積分奇異點(diǎn)處理［9，12］、邊界場連續(xù)性問題［12］等開展了大量的研究。一般來說，解析法表達(dá)式比較復(fù)雜、適應(yīng)性較差、計(jì)算量較大。數(shù)值法主要通過有限差分法［13］、有限體積法［14］、有限單元法［15-16］求解重磁異常滿足的泊松方程，但是當(dāng)計(jì)算大量位置點(diǎn)的異常場時(shí)，耗時(shí)較長。

頻率域方法的原理是采用數(shù)值方法對場源在某條測線（或平面／三維空間）產(chǎn)生的重磁異常頻譜（傅里葉變換解析表達(dá)式）進(jìn)行計(jì)算，并通過反傅里葉變換得到異常場［17］。其優(yōu)點(diǎn)是表達(dá)式簡單，計(jì)算效率相對較高，但是傅里葉變換的截?cái)嘈?yīng)限制了該方法的普遍適用性。Bhattacharyya［3］推導(dǎo)了任意磁化方向的長方體磁場頻譜表達(dá)式；Pedersen［18］給出了直立圓柱體和多面體模型的重磁異常頻率域表達(dá)式；熊光楚［19］給出了長方體模型重力位的三維傅里葉變換式；后來一些學(xué)者對物性連續(xù)變化模型的頻率域 解 析 式［20-24］、Parker模 型 頻 率 域 表 達(dá) 式［25-26］、偏移抽樣方法提高反傅里葉變換數(shù)值精度［27］開展了大量研究。Tontini等［28］實(shí)現(xiàn)了重磁異常三維全空間頻率域正演，并研究了FFT 擴(kuò)邊法與誤差的關(guān)系；Wu等［29-30］利用高斯FFT 提高了反傅里葉變換的數(shù)值精度，降低了由于FFT 引起的強(qiáng)制周期化、邊界震蕩效應(yīng)等影響，并研究了地形對重磁異常的影響，分析了物性連續(xù)變化模型的梯度場；Ren等［31］提出一種基于非結(jié)構(gòu)化網(wǎng)格的自適應(yīng)快速多極子重磁場及其梯度張量正演方法，借助并行計(jì)算實(shí)現(xiàn)了快速、高精度的三維正演。

目前，隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展和重磁異常正演算法的不斷改進(jìn)，頻率域方法以其簡捷、準(zhǔn)確和高效等優(yōu)勢在數(shù)值模擬中的應(yīng)用越來越廣泛。但是，對于面向?qū)嶋H應(yīng)用的超大規(guī)模復(fù)雜三維模型數(shù)值模擬問題，由于網(wǎng)格剖分單元數(shù)目劇增，計(jì)算量和存儲量巨大，常規(guī)的空間域和頻率域重磁異常正演方法已很難同時(shí)滿足精度與效率上的要求。為此，本文提出一種重磁位場三維數(shù)值模擬方法，充分利用一維形函數(shù)積分的高精度性、不同波數(shù)之間高度并行性及傅里葉變換的高效性，實(shí)現(xiàn)高效、高精度、大規(guī)模的重磁位場數(shù)值模擬。

1 理論與方法

弱磁性體的磁位V 及重力位φ 的空間域積分表達(dá)分別為［32］

式中：M 為空間域磁化強(qiáng)度矢量；ρ為剩余密度；μ為真空磁導(dǎo)率，取值4π×10－7H／m；D 為萬有引力常量，取值6.674×10－11m3kg－1s－2；（x，y，z）為觀測點(diǎn)坐標(biāo)；（x′，y′，z′）為異常點(diǎn)坐標(biāo)；為哈密頓算符；G 為格林函數(shù)。

對式（1）沿x、y 方向做二維傅里葉變換，分別得到磁力位和重力位一維積分公式

式（2）的推導(dǎo)需利用波數(shù)域格林函數(shù)，其表達(dá)式為

磁位和重力位的一階導(dǎo)數(shù)分別為

對應(yīng)的二階導(dǎo)數(shù)分別為

式中下標(biāo)“B”和“g”分別代表磁力場和重力場。

利用傅里葉變換的微分性質(zhì)，對式（4）和式（5）做二維傅里葉變換，得到波數(shù)域重磁異常場和張量的表達(dá)式分別為

式（2）為重磁位垂向一維積分表達(dá)式。從該表達(dá)式可以看出，不同波數(shù)一維空間域積分是相互獨(dú)立的，可以并行計(jì)算；該一維積分深度方向可離散為多個(gè)單元積分之和，每個(gè)積分單元采用形函數(shù)法插值計(jì)算，得到單元積分的解析表達(dá)式，可提高計(jì)算精度和效率。

2 基于形函數(shù)的一維積分

本文采用基于二次插值的形函數(shù)法［31］計(jì)算傅里葉變換后的一維積分。將一維積分沿垂直方向剖分為m 個(gè)單元，采用二次形函數(shù)描述物性參數(shù)變化。從式（2）可以看出，重力位與磁位的積分表達(dá)式相似，只是系數(shù)不同；重力異常與磁異常積分中，不同方向的積分表達(dá)式也相似，只是系數(shù)不同。因此，下面以x 方向磁位積分為例，介紹形函數(shù)法計(jì)算一維積分的具體過程。

對積分區(qū)域進(jìn)行離散

式中i表示積分單元。

對磁化強(qiáng)度在節(jié)點(diǎn)處采用二次形函數(shù)插值

式中：N1、N2、N3為二次形函數(shù)；i1、i2、i3分別為積分單元i的3個(gè)插值節(jié)點(diǎn)。

將式（9）代入式（8）中，整理得

由式（10）離散后的一維單元積分可推導(dǎo)出解析表達(dá)式，該表達(dá)式計(jì)算精度較高，接近于解析解。對于物性參數(shù)滿足二次變化的復(fù)雜形體，垂向單元剖分間隔可根據(jù)需要靈活選擇。

3 算法流程

本文按照如下流程進(jìn)行重磁位場數(shù)值模擬計(jì)算。

（1）網(wǎng)格剖分。給定計(jì)算區(qū)域，確定沿x、y、z方向的剖分單元個(gè)數(shù)C1、C2、C3。水平方向均勻剖分，垂直方向可任意剖分，得到三維離散坐標(biāo)點(diǎn)（xm，yn，zl），給定計(jì)算的觀測平面z。

（2）模型參數(shù)設(shè)定。給所有的剖分單元賦予剩余密度或磁化率值。

（3）離散波數(shù)計(jì)算。由FFT 法的波數(shù)計(jì)算公式得出水平離散坐標(biāo)（xm，yn）的離散波數(shù)（kx，ky）。

（4）重磁空間域二維傅里葉變換。對空間域重磁位場剩余密度或者磁化強(qiáng)度進(jìn)行水平二維離散傅里葉變換。

（5）積分計(jì)算。采用形函數(shù)法計(jì)算垂向一維積分，得到波數(shù)譜。

（6）重磁波數(shù)譜二維傅里葉反變換。利用FFT算法，對波數(shù)譜場或張量進(jìn)行二維離散傅里葉反變換，得到空間域重磁位場或張量。

4 算例分析

4.1 正確性驗(yàn)證

設(shè)計(jì)一個(gè)棱柱體模型，采用基于四個(gè)高斯點(diǎn)的高斯—FFT法［28］計(jì)算重磁異常場及張量的數(shù)值解，并與模型解析解［32］進(jìn)行對比，驗(yàn)證該算法的正確性。

模型如圖1所示。異常體大小為200m×200m×200m，頂面埋深h為200m。觀測面為z＝0的水平面，模擬區(qū)域?yàn)?00m×800m×400m，觀測點(diǎn)個(gè)數(shù)為201×201。棱柱異常體與模擬區(qū)域水平中心點(diǎn)重合。模擬區(qū)域網(wǎng)格節(jié)點(diǎn)個(gè)數(shù)為201×201×201，空間采樣間隔均為4m。設(shè)定異常體的磁化強(qiáng)度為0.01A／m，背景磁場為4500n T，磁傾角為45°，磁偏角為5°；設(shè)定其剩余密度為1000kg／m3。

圖1 棱柱體模型示意圖

圖2為磁異常場在z＝0平面的數(shù)值解、解析解與絕對誤差?？梢钥闯?，磁異常場的絕對誤差最大約為5×10－3n T。圖3為磁異常梯度張量在z＝0平面的數(shù)值解、解析解與絕對誤差，可以看出，磁異常梯度張量的絕對誤差最大約為5×10－5n T／m。磁異常場和磁異常梯度張量的絕對誤差均比相應(yīng)解析解小三個(gè)數(shù)量級。

圖4為重力異常場在z＝0平面的數(shù)值解、解析解與絕對誤差?？梢钥闯?，重力異常場的絕對誤差最大約為2×10－4mGal。圖5為重力異常梯度張量在z＝0平面的數(shù)值解、解析解與絕對誤差，可以看出，重力異常梯度張量的絕對誤差最大約為0.02E。同樣地，重力異常場和重力異常梯度張量的絕對誤差值也均比相應(yīng)解析解低三個(gè)數(shù)量級。重磁場數(shù)值模擬結(jié)果表明，本文算法正確，且精度高。

圖2 z＝0平面磁異常場的數(shù)值解（左）、解析解（中）及絕對誤差（右）

4.2 復(fù)雜形體數(shù)值模擬適用性研究

對波數(shù)域的一維積分采用二次形函數(shù)法可得到每個(gè)單元的近似解析解，該方法計(jì)算精度較高，適用于復(fù)雜介質(zhì)的數(shù)值模擬計(jì)算。以磁異常數(shù)值模擬為例設(shè)計(jì)一個(gè)復(fù)雜形體模型。異常體的水平方向磁化率不變，垂直方向磁化率呈二次變化。垂直方向分別采用均勻剖分與形函數(shù)插值兩種方法獲得數(shù)值解。通過數(shù)值解與解析解［25］的對比分析，可驗(yàn)證本文方法對于復(fù)雜模型的適用性。為了研究統(tǒng)計(jì)整個(gè)觀測面的誤差情況，引入相對均方根誤差［30］

圖3 z＝0平面磁梯度張量的數(shù)值解（左）、解析解（中）及絕對誤差（右）

圖4 z＝0平面重力異常場的數(shù)值解（左）、解析解（中）及絕對誤差（右）

式中：P 和N 分別是x 和y 方向的節(jié)點(diǎn)數(shù)；Bij是磁場（張量）數(shù)值解；^Bij是磁場（張量）解析解。該誤差統(tǒng)計(jì)方式能突出大異常值所占的比重。設(shè)計(jì)兩個(gè)不同模型，水平方向磁化率均不變，垂向磁化率從0遞增到0.02SI，其中模型1（圖6）垂直方向磁化率變化為0.02／60002（z－2000）2（2000≤z≤8000），模型2（圖7）磁化率變化為0.02／20002（z－4000）2（4000≤z≤6000），水平方向采樣間距均為100m。兩個(gè)模型的背景磁感應(yīng)強(qiáng)度均為45000n T，磁傾角為45°，磁偏角為5°。

圖8和圖9分別為模型1磁場分量和磁張量隨著采樣間距變化的數(shù)值解與解析解的相對均方根誤差變化圖?？梢钥闯觯瑑煞N方法在垂直采樣間距為25m 時(shí)的計(jì)算精度相近，在積分單元內(nèi)磁化率均勻剖分的誤差隨著采樣間距的增大而增大，而形函數(shù)二次插值的計(jì)算精度基本保持不變。

圖10和圖11分別為模型2磁場分量和磁張量隨著采樣間距變化的數(shù)值解與解析解的相對均方根誤差變化圖。

與模型算例1計(jì)算結(jié)果對比可以看出，兩個(gè)模型在兩種剖分方法下的誤差變化趨勢基本一致。但是與模型1相比，模型2中單元內(nèi)磁化率均勻剖分的計(jì)算精度有所降低，這是由于模型2中模型介質(zhì)磁化率的變化更劇烈，而單元內(nèi)磁化率形函數(shù)二次插值的計(jì)算精度變化較小。這再次驗(yàn)證了該方法對于復(fù)雜介質(zhì)模型具有較好的適用性。

對于模型區(qū)域剖分網(wǎng)格單元數(shù)為100×100×100、觀測點(diǎn)數(shù)為101×101的三維正演計(jì)算，傳統(tǒng)空間域累加算法耗時(shí)約800s，波數(shù)域Parker公式計(jì)算結(jié)果約耗時(shí)8.4s，本文算法約耗時(shí)5.1s，加速比分別為156.8與1.65，驗(yàn)證了本算法的高效性。

圖5 z＝0平面重力梯度張量的數(shù)值解（左）、解析解（中）及絕對誤差（右）

圖6 模型1垂向磁化率變化趨勢圖

圖7 模型2垂向磁化率變化趨勢圖

圖8 模型1磁場分量相對均方根誤差統(tǒng)計(jì)曲線

圖9 模型1磁張量相對均方根誤差統(tǒng)計(jì)曲線

圖10 模型2磁場分量相對均方根誤差

圖11 模型2磁張量相對均方根誤差

5 結(jié)論

本文提出一種高效、高精度的重磁位場三維數(shù)值模擬方法。該方法沿水平方向進(jìn)行二維傅里葉變換，把重磁位場三維積分轉(zhuǎn)化為垂向一維積分問題。對于離散后的單元積分采用形函數(shù)二次插值計(jì)算，可得出單元積分的解析表達(dá)式，理論上具有較高的計(jì)算精度。

設(shè)計(jì)棱柱體模型，模型的解析解與本文方法的數(shù)值解對比表明，本文提出的重磁位場數(shù)值模擬方法正確且計(jì)算精度高。

設(shè)計(jì)復(fù)雜連續(xù)變化介質(zhì)模型，模型解析解與數(shù)值解對比表明，本文提出的形函數(shù)插值計(jì)算一維積分與傳統(tǒng)的棱柱體均勻剖分相比具有更高的計(jì)算精度，更適用于變化復(fù)雜的實(shí)際介質(zhì)的模擬計(jì)算。

附錄A 形函數(shù)單元積分系數(shù)推導(dǎo)過程

對第i個(gè)單元進(jìn)行積分時(shí)，根據(jù)觀測點(diǎn)與積分單元的相對位置，分三種情況進(jìn)行討論：

（1）當(dāng)觀測點(diǎn)在該單元上方時(shí)，即z＜z′，式（A-1）為

（2）當(dāng)觀測點(diǎn)在該單元下方時(shí)，即z＞z′i＋2，式（A-1）為

（3）當(dāng)觀測點(diǎn)在該單元內(nèi)部時(shí)，即z′≤z≤z′i＋2，式（A-1）為

令γ＝k，式（A-2）變?yōu)?/p>

其中

得到單元積分

求解過程如下式（A-9）的單元形函數(shù)展開與式（A-6）類似，只是積分的上、下限不同。

利用每個(gè)單元的插值函數(shù)，可將所求積分表示為