復雜形體重磁位場三維高效高精度數值模擬

2020-10-17 07:44:30周印明戴世坤凌嘉宣胡曉穎

石油地球物理勘探 2020年5期

關鍵詞：模型

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

（①中南大學地球科學與信息物理學院，湖南長沙410083；②中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室，湖南長沙410083；③東方地球物理公司綜合物化探處，河北涿州072751；④桂林理工大學地球科學學院，廣西桂林541006）

0 引言

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

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

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

1 理論與方法

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

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

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

式（2）的推導需利用波數域格林函數，其表達式為

磁位和重力位的一階導數分別為

對應的二階導數分別為

式中下標“B”和“g”分別代表磁力場和重力場。

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

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

2 基于形函數的一維積分

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

對積分區域進行離散

式中i表示積分單元。

對磁化強度在節點處采用二次形函數插值

式中：N1、N2、N3為二次形函數；i1、i2、i3分別為積分單元i的3個插值節點。

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

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

3 算法流程

本文按照如下流程進行重磁位場數值模擬計算。

（1）網格剖分。給定計算區域，確定沿x、y、z方向的剖分單元個數C1、C2、C3。水平方向均勻剖分，垂直方向可任意剖分，得到三維離散坐標點（xm，yn，zl），給定計算的觀測平面z。

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

（3）離散波數計算。由FFT 法的波數計算公式得出水平離散坐標（xm，yn）的離散波數（kx，ky）。

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

（5）積分計算。采用形函數法計算垂向一維積分，得到波數譜。

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

4 算例分析

4.1 正確性驗證

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

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

圖1 棱柱體模型示意圖

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

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

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

4.2 復雜形體數值模擬適用性研究

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

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

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

式中：P 和N 分別是x 和y 方向的節點數；Bij是磁場（張量）數值解；^Bij是磁場（張量）解析解。該誤差統計方式能突出大異常值所占的比重。設計兩個不同模型，水平方向磁化率均不變，垂向磁化率從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。兩個模型的背景磁感應強度均為45000n T，磁傾角為45°，磁偏角為5°。

圖8和圖9分別為模型1磁場分量和磁張量隨著采樣間距變化的數值解與解析解的相對均方根誤差變化圖。可以看出，兩種方法在垂直采樣間距為25m 時的計算精度相近，在積分單元內磁化率均勻剖分的誤差隨著采樣間距的增大而增大，而形函數二次插值的計算精度基本保持不變。

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

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

對于模型區域剖分網格單元數為100×100×100、觀測點數為101×101的三維正演計算，傳統空間域累加算法耗時約800s，波數域Parker公式計算結果約耗時8.4s，本文算法約耗時5.1s，加速比分別為156.8與1.65，驗證了本算法的高效性。