胡正旺 杜勁松*③ 孫石達(dá) 陳 超

(①中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北武漢 430074；②地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430074；③地質(zhì)過(guò)程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430074；④河北地質(zhì)大學(xué)勘查技術(shù)與工程學(xué)院，河北石家莊 050031)

0 引言

在地球物理學(xué)方法體系中，磁法發(fā)展最早，至今仍被廣泛應(yīng)用于地質(zhì)調(diào)查、資源勘探、軍事與國(guó)防建設(shè)、工程勘察與環(huán)境監(jiān)測(cè)等諸多領(lǐng)域[1-3]。地磁場(chǎng)測(cè)量的平臺(tái)、方式以及物理量多種多樣，目前磁法勘探普遍采用標(biāo)量磁力儀，基于航空、地面與船載等移動(dòng)平臺(tái)測(cè)量地磁場(chǎng)總強(qiáng)度的空間分布[4]。隨著地磁場(chǎng)測(cè)量精度和測(cè)點(diǎn)定位精度的大幅提高，以及計(jì)算機(jī)與計(jì)算技術(shù)的快速發(fā)展，對(duì)深部目標(biāo)及小型磁性體的探測(cè)需求日益增長(zhǎng)。以往在磁力異常的校正計(jì)算、數(shù)據(jù)預(yù)處理、正演、反演及定量解釋中的許多假設(shè)條件和近似條件逐漸變得不再合理，在傳統(tǒng)算法中被忽略的系列因素也逐步納入考慮范圍，例如自退磁效應(yīng)[5]、磁化率各向異性[6]、剩余磁化強(qiáng)度[7]、總磁場(chǎng)強(qiáng)度異常的非線性效應(yīng)[8-17]等。

總磁場(chǎng)強(qiáng)度異常(ΔT)的非線性效應(yīng)，指的是實(shí)際測(cè)量得到的磁異常ΔT為地磁場(chǎng)總場(chǎng)(T)的模量|T|與主磁場(chǎng)(T0)的模量|T0|之差(即ΔT=|T|-|T0|)，與磁場(chǎng)矢量異常(Ta)以及巖石圈磁源的磁矩不滿足線性關(guān)系。在實(shí)際數(shù)據(jù)預(yù)處理、正反演和定量解釋時(shí)，通常將Ta在主磁場(chǎng)或正常地磁場(chǎng)T0方向上的投影分量ΔTPro近似為ΔT[10]，這就意味著ΔT與其他磁分量異常一樣，與場(chǎng)源的磁矩呈簡(jiǎn)單的線性關(guān)系，與主磁場(chǎng)T0無(wú)耦合關(guān)系，并且在場(chǎng)源外部空間滿足拉普拉斯方程從而具備調(diào)和性質(zhì)，因而可以簡(jiǎn)化實(shí)測(cè)ΔT數(shù)據(jù)的預(yù)處理、正反演及定量解釋過(guò)程[8-9]。該近似的前提是|T0|?|Ta|，例如：當(dāng)|T0|=30000nT、|Ta|=500nT時(shí)，ΔT近似為ΔTPro的誤差僅為4nT；當(dāng)|T0|=70000nT、|Ta|=10000nT時(shí)，ΔT近似為ΔTPro的誤差可達(dá)714nT[11]；隨著磁異常模量Ta=|Ta|的增大，誤差按照近似指數(shù)關(guān)系迅速增大[8-9]。總之，當(dāng)Ta幅值較大時(shí)，這種近似處理會(huì)產(chǎn)生較大的誤差，進(jìn)而對(duì)實(shí)測(cè)ΔT的處理(例如空間延拓曲化平、分量與其他參量的轉(zhuǎn)換、化極、磁源重力異常、導(dǎo)數(shù)計(jì)算等)、正反演及定量解釋產(chǎn)生較大的影響。

針對(duì)該問(wèn)題，Zhen等[12]和甄慧翔等[13-14]在波數(shù)域提出了空間域基于ΔT精確計(jì)算ΔTPro的最優(yōu)化方法，Sun等[15]和孫石達(dá)等[16]提出基于三維成像反演的非線性等效源方法。這兩種方法均需要求解大型方程組，不僅所需內(nèi)存較大且計(jì)算效率較低，限制了其在實(shí)際生產(chǎn)與研究中的廣泛應(yīng)用。針對(duì)此問(wèn)題，本文提出一種高效率的迭代方法，即采用ΔT近似為ΔTPro的理論誤差公式進(jìn)行迭代計(jì)算，并通過(guò)模型試驗(yàn)和實(shí)際數(shù)據(jù)應(yīng)用驗(yàn)證了該方法的準(zhǔn)確性、穩(wěn)定性與計(jì)算效率，旨在推動(dòng)高精度實(shí)測(cè)ΔT數(shù)據(jù)的實(shí)用化進(jìn)程。

1 方法原理與算法流程

首先介紹由ΔT迭代計(jì)算ΔTPro的方法原理，然后給出算法流程，最后從理論上分析該迭代算法計(jì)算精度與計(jì)算效率的影響因素。

圖1所示為總磁場(chǎng)強(qiáng)度異常ΔT、磁異常模量Ta與磁場(chǎng)矢量異常Ta在主磁場(chǎng)T0方向投影ΔTPro之間的關(guān)系示意圖。若忽略外源變化場(chǎng)及其在地球內(nèi)部感應(yīng)生成的二次磁場(chǎng)，地磁場(chǎng)T可視為主磁場(chǎng)T0與磁場(chǎng)矢量異常Ta之矢量和

T=T0+Ta

(1)

如引言所述，實(shí)測(cè)總磁場(chǎng)強(qiáng)度異常ΔT往往表示為總磁場(chǎng)強(qiáng)度|T|與主磁場(chǎng)強(qiáng)度|T0|(即T0)之差

ΔT=|T|-|T0|

(2)

而在實(shí)際數(shù)據(jù)預(yù)處理、正反演及定量解釋時(shí)，為了簡(jiǎn)化計(jì)算，通常將ΔTPro近似為ΔT(圖1)，即

圖1 總磁場(chǎng)強(qiáng)度異常(ΔT)、磁異常模量(Ta)與磁場(chǎng)矢量異 常Ta在主磁場(chǎng)T0方向投影ΔTPro之間的關(guān)系示意圖

ΔTPro=|Ta|cosθ

(3)

顯然，將ΔTPro近似表示為ΔT存在誤差

E=ΔT-ΔTPro

(4)

根據(jù)圖1，經(jīng)過(guò)推導(dǎo)與化簡(jiǎn)可得誤差E的表達(dá)式[8]為

(5)

由前文所述的假設(shè)條件Ta≥ΔT，可得E≥0、 ΔT≥ΔTPro[8]。

聯(lián)合式(4)與式(5)可得

(6)

式(6)即為基于ΔT精確計(jì)算投影分量ΔTPro的方法基礎(chǔ)。若假設(shè)研究區(qū)主磁場(chǎng)T0的大小和方向不變，根據(jù)投影分量ΔTPro可在波數(shù)域快速轉(zhuǎn)換計(jì)算磁場(chǎng)矢量異常Ta的模Ta，又因ΔT為實(shí)測(cè)數(shù)據(jù)，則根據(jù)式(6)可以構(gòu)建方程組，進(jìn)而采用最優(yōu)化方法求解ΔTPro。這是甄慧翔等[12-14]提出的方法思路，該方法需求解大型方程組，因而計(jì)算效率較低。

本文基于式(6)，提出采用迭代算法求解ΔTPro，即

(7)

式中：i為非負(fù)整數(shù)，表示迭代次數(shù)；主磁場(chǎng)強(qiáng)度T0在研究區(qū)內(nèi)為常數(shù)。ΔT(i)的計(jì)算公式為

(8)

(9)

(10)

式中：F與F-1分別表示二維傅里葉正變換與反變換；j為虛數(shù)單位；L0、M0與N0分別為主磁場(chǎng)單位矢量t0與x、y、z軸之間的方向余弦；u、v分別表示x、y方向的波數(shù)。因此，根據(jù)式(9)可得

(11)

ε=max(|ΔT-ΔT(i)|)≤ε0

(12)

圖2 本文迭代算法流程圖 圖中虛線表示迭代更新過(guò)程

(1)對(duì)實(shí)測(cè)數(shù)據(jù)ΔT進(jìn)行數(shù)據(jù)預(yù)處理，依次進(jìn)行奇異值剔除、數(shù)據(jù)網(wǎng)格化、擴(kuò)邊、初步去噪；

(2)給定ε0；

由上述方法原理與算法流程可以看出，與Sun等[15]和孫石達(dá)等[16]采用的基于三維成像反演的非線性等效源方法不同，本文方法與文獻(xiàn)[12-14]所提出的方法在原理上是一致的，差異僅在于所采用的計(jì)算方法不同：本文采用了迭代算法，后者采用的是基于最優(yōu)化方法求解線性方程組。甄慧翔[13]明確指出：基于式(7)計(jì)算ΔTPro的精度影響因素僅在于采用傅里葉變換計(jì)算磁場(chǎng)分量異常的過(guò)程，而與場(chǎng)源相關(guān)的剩余磁化效應(yīng)、自退磁效應(yīng)、磁性結(jié)構(gòu)復(fù)雜性等均無(wú)關(guān)系；對(duì)于實(shí)測(cè)ΔT數(shù)據(jù)所含誤差的影響，由式(9)可知，分量轉(zhuǎn)換的波數(shù)域轉(zhuǎn)換因子并非如空間延拓和導(dǎo)數(shù)計(jì)算那樣對(duì)噪聲具有明顯的壓制或放大作用，而是略微受主磁場(chǎng)方向的影響，甄慧翔[13]所做的模型試驗(yàn)、以及Coleman等[18]基于等效源的相關(guān)研究也證明了此點(diǎn)。此外，從上述迭代計(jì)算過(guò)程還可以看出，計(jì)算速度主要決定于迭代總次數(shù)以及每次迭代所需的1次傅里葉正變換和3次傅里葉反變換。因此，本文所提方法的計(jì)算精度與計(jì)算效率在根源上主要取決于式(9)所示的傅里葉正、反變換。

2 模型試驗(yàn)

為了檢驗(yàn)本文方法的可靠性與穩(wěn)定性，本節(jié)基于理論磁性結(jié)構(gòu)模型、通過(guò)正演計(jì)算模擬觀測(cè)數(shù)據(jù)與理論計(jì)算結(jié)果，對(duì)模擬的觀測(cè)數(shù)據(jù)采用本文方法進(jìn)行計(jì)算，并對(duì)計(jì)算結(jié)果與理論結(jié)果進(jìn)行對(duì)比和分析。

2.1 理論模型

鑒于本文方法與甄慧翔[13]所用方法的原理相似，為了避免重復(fù)性工作，不再進(jìn)行數(shù)據(jù)噪聲與背景場(chǎng)、剩余磁化效應(yīng)、自退磁效應(yīng)、磁性結(jié)構(gòu)復(fù)雜性等模型試驗(yàn)研究，相關(guān)內(nèi)容可參照文獻(xiàn)[12-14]。本文重點(diǎn)研究三個(gè)方面的內(nèi)容：一是分析迭代是否能夠收斂，若能夠收斂，進(jìn)一步分析其計(jì)算精度；二是分析計(jì)算精度與迭代收斂速度的影響因素；三是分析計(jì)算效率。

本文針對(duì)性地設(shè)計(jì)了四組模型試驗(yàn)，均采用簡(jiǎn)單的磁偶極源，場(chǎng)源參數(shù)、主磁場(chǎng)參數(shù)與觀測(cè)數(shù)據(jù)參數(shù)見(jiàn)表1，所有試驗(yàn)均在直角坐標(biāo)系中進(jìn)行，x軸指向北、y軸指向東、z軸指向地下，觀測(cè)面均為z=0的水平面，且僅考慮感應(yīng)磁化問(wèn)題，磁偶極源磁場(chǎng)的計(jì)算公式參見(jiàn)文獻(xiàn)[19]中的式(32)～式(34)。第1組模型試驗(yàn)用于測(cè)試方法的有效性，即迭代是否收斂并分析計(jì)算精度；第2組和第3組模型試驗(yàn)分別改變磁偶極源的磁矩及主磁場(chǎng)的強(qiáng)度、傾角和偏角，分析磁場(chǎng)矢量異常的模量大小、主磁場(chǎng)強(qiáng)度和方向?qū)κ諗啃院陀?jì)算精度的影響；第4組模型試驗(yàn)的正演觀測(cè)數(shù)據(jù)采用不同的點(diǎn)距與線距，主要測(cè)試分析計(jì)算效率隨數(shù)據(jù)總量的變化規(guī)律。

表1 理論磁性模型和磁場(chǎng)觀測(cè)參數(shù)表

2.2 試驗(yàn)結(jié)果及其分析

2.2.1 方法有效性

第1組模型試驗(yàn)結(jié)果見(jiàn)圖3。由圖可見(jiàn)，擬合差ε隨迭代次數(shù)增加呈先快速衰減、然后趨于穩(wěn)定的特征，迭代2次時(shí)最大擬合誤差即可達(dá)到約10nT，可見(jiàn)本文迭代方法具有快速與穩(wěn)定的收斂特征。此外，圖3中ε趨于常值，說(shuō)明迭代終止條件采用式(12)定義的閾值ε0不再合適，而應(yīng)采用相鄰兩次迭代的擬合誤差之差δ=|εi+1-εi|更合理。若ε0設(shè)置太小，則可能致使迭代結(jié)果永遠(yuǎn)達(dá)不到誤差要求，因此在實(shí)際數(shù)據(jù)處理時(shí)可以采用δ作為迭代終止條件，也可以通過(guò)設(shè)定最大迭代次數(shù)停止迭代計(jì)算過(guò)程。

Zhen等[12]和甄慧翔等[13-14]利用式(7)進(jìn)行迭代時(shí)，ΔT(i)采用了實(shí)測(cè)ΔT。從圖3可見(jiàn)，ΔT(i)采用式(8)計(jì)算結(jié)果與采用實(shí)測(cè)ΔT兩種策略對(duì)迭代收斂性的影響沒(méi)有顯著差異。但是筆者認(rèn)為，式(5)所示誤差E與ΔT、T0、Ta之間是耦合的， 應(yīng)采用式(8)計(jì)算，若采用實(shí)測(cè)ΔT數(shù)據(jù)代替可能會(huì)存在迭代過(guò)程振蕩的風(fēng)險(xiǎn)，盡管在此例中并未出現(xiàn)這種現(xiàn)象。

圖3 第1組模型試驗(yàn)擬合誤差ε隨迭代次數(shù)的變化曲線

圖4為第1組模型的理論磁異常分布，圖5為理論值與本文方法計(jì)算結(jié)果的差值，其統(tǒng)計(jì)參數(shù)見(jiàn)表2?？梢?jiàn)轉(zhuǎn)換處理的計(jì)算誤差整體較低，其中磁場(chǎng)矢量異常的模量Ta計(jì)算誤差相對(duì)較大，但也僅約±7nT，相比其最大幅值10614nT，該誤差完全可以忽略，因此本文算法的計(jì)算精度較高。若直接將ΔT近似為ΔTPro，采用式(9)與式(11)分別計(jì)算三分量異常和模量異常，其計(jì)算誤差見(jiàn)圖6和表2，其中Za分量的最大誤差達(dá)到了958.9nT，可見(jiàn)傳統(tǒng)方法在高磁環(huán)境下具有較大的轉(zhuǎn)換計(jì)算誤差。

圖4 第1組模型的理論磁異常圖 (a)ΔT； (b)ΔTPro； (c)Ha,x； (d)Ha,y； (e)Za； (f)Ta 等值線間隔為1000nT；黑色虛線、點(diǎn)劃線與實(shí)線分別表示負(fù)值、零值與正值。圖5、圖6、圖9～圖11、圖14、圖15同

圖5 第1組模型的理論值與本文方法計(jì)算結(jié)果的差值 (a)ΔT； (b)ΔTPro； (c)Ha,x； (d)Ha,y； (e)Za； (f)Ta 等值線間隔為0.5nT

圖3中另一個(gè)明顯特征是當(dāng)?shù)_(dá)到一定次數(shù)(此例為3次)后，擬合誤差ε趨于穩(wěn)定(此例約趨于5nT)，其原因可以從如下兩個(gè)方面進(jìn)行分析。

圖6 第1組模型的理論值與傳統(tǒng)方法計(jì)算值的差值 (a)Ha,x； (b)Ha,y； (c) Za； (d)Ta 等值線間隔為50nT

表2 第1組模型試驗(yàn)的理論值、計(jì)算值及誤差統(tǒng)計(jì)表

(1)有限項(xiàng)傅里葉級(jí)數(shù)展開(kāi)本身的Gibbs現(xiàn)象可能會(huì)產(chǎn)生邊界效應(yīng)。離散傅里葉變換的前提是信號(hào)在計(jì)算區(qū)域內(nèi)是周期性的，顯然此例所示磁力異常數(shù)據(jù)在研究區(qū)域內(nèi)不具備完全的周期性，因?yàn)榇帕Ξ惓Ｔ谶吔绮荒芡耆嗟?，因而去除常值之后在邊界處存在不連續(xù)點(diǎn)。如此采用周期函數(shù)擬合非周期信號(hào)，必然在邊界區(qū)域產(chǎn)生振蕩，其壓制方法是擴(kuò)大數(shù)據(jù)范圍[20-25]。一方面在計(jì)算完成之后可通過(guò)截取目標(biāo)區(qū)域結(jié)果避免邊界振蕩；另一方面可通過(guò)合理擴(kuò)邊，使得邊界處的數(shù)值趨于相等，從而壓制邊界效應(yīng)。事實(shí)上，在此組模型試驗(yàn)中，已經(jīng)先將數(shù)據(jù)向南北和東西方向各擴(kuò)邊了50m，邊界處的數(shù)值已經(jīng)均接近于零值，所以邊界效應(yīng)并不明顯。但是，擬合誤差隨著迭代次數(shù)增加依然不能趨于零值，這說(shuō)明除了邊界效應(yīng)之外還存在其他原因。

(2)觀察表2可以發(fā)現(xiàn)，計(jì)算誤差實(shí)際上存在常值成分，這源于有限項(xiàng)傅里葉級(jí)數(shù)展開(kāi)的常值項(xiàng)，舒晴等[26]稱之為波數(shù)域奇點(diǎn)，認(rèn)為這是轉(zhuǎn)換誤差產(chǎn)生的根源，可通過(guò)減小點(diǎn)距、擴(kuò)大計(jì)算范圍有效提高轉(zhuǎn)換處理的精度，其中擴(kuò)大計(jì)算范圍的改善效果非常顯著。

甄慧翔等[14]所做的疊加背景場(chǎng)模型試驗(yàn)也已證明，背景場(chǎng)的存在使基于快速傅里葉變換的分量轉(zhuǎn)換計(jì)算產(chǎn)生較大的誤差，其本質(zhì)原因即在于上述的邊界效應(yīng)和波數(shù)域奇點(diǎn)問(wèn)題。因此，在實(shí)際應(yīng)用中，建議盡量采用較大的計(jì)算區(qū)域并事先適當(dāng)消除區(qū)域均值或者趨勢(shì)場(chǎng)，當(dāng)然這是在一般的情況下，即關(guān)注的目標(biāo)是局部異常。另外，波數(shù)域奇點(diǎn)問(wèn)題可以采用偏移抽樣理論[27-30]進(jìn)行壓制，這是今后待改進(jìn)之所在，本文不再過(guò)多論及。

2.2.2 迭代收斂速度與計(jì)算精度的影響因素

由式(5)可知，采用ΔT近似ΔTPro所產(chǎn)生的誤差主要源于Ta和T0，因此設(shè)計(jì)表1所示的第2組與第3組模型試驗(yàn)，分別調(diào)查二者對(duì)轉(zhuǎn)換計(jì)算精度和迭代收斂速度的影響。

由于影響Ta大小的主要因素包括磁偶極源的磁矩大小及場(chǎng)源埋深，而埋深變化會(huì)影響磁異常在測(cè)區(qū)的分布形態(tài)，因此第2組模型試驗(yàn)僅考慮磁矩大小的影響。兩個(gè)模型對(duì)應(yīng)的擬合誤差ε隨迭代次數(shù)的變化曲線見(jiàn)圖7?？梢钥闯?，場(chǎng)源磁矩越大，迭代收斂速度越慢、計(jì)算精度越低。若取δ=1nT，則兩個(gè)模型達(dá)到該閾值所需的最小迭代次數(shù)分別為4與7。分析其原因在于：磁場(chǎng)矢量異常的模量越大，式(7)的逼近誤差也越大，因此迭代收斂變慢；并且在擴(kuò)邊相同的情況下，研究區(qū)域邊界處的數(shù)值越高，波數(shù)域奇點(diǎn)的影響也就更加嚴(yán)重，因而計(jì)算精度越低。

圖8為第3組模型的擬合誤差ε隨迭代次數(shù)的變化曲線。需要說(shuō)明的是，此組模型的主磁場(chǎng)強(qiáng)度和方向以及磁偶極源的磁化方向均是變化的，但是通過(guò)改變磁偶極源體積保持了磁偶極源的磁矩不變，即1×106A·m2。由圖8可見(jiàn)，主磁場(chǎng)的強(qiáng)度T0和方向(I與D)對(duì)迭代收斂速度的影響均較小(圖中三條曲線在第1次迭代后近于平行)，差異僅在于計(jì)算精度，即當(dāng)磁化傾角與主磁場(chǎng)傾角為90°時(shí)，無(wú)論主磁場(chǎng)強(qiáng)度為多大，其計(jì)算誤差均較大且均趨于相同的擬合誤差ε；而當(dāng)磁傾角與主磁傾角接近0°時(shí)，則具有較高的計(jì)算精度。分析其原因在于：主磁場(chǎng)強(qiáng)度T0對(duì)式(7)的逼近誤差的影響較小，因而主磁場(chǎng)強(qiáng)度T0對(duì)于迭代收斂性的影響程度較低；對(duì)于計(jì)算精度，當(dāng)磁傾角與主磁場(chǎng)傾角為90°時(shí)，在擴(kuò)邊相同的情況下，研究區(qū)域內(nèi)ΔT的均值越高，波數(shù)域奇點(diǎn)的影響也就越嚴(yán)重。

第2組與第3組模型試驗(yàn)轉(zhuǎn)換計(jì)算結(jié)果及其與理論值的差值分別見(jiàn)圖9和圖10。由圖可見(jiàn)，第2組模型的ΔTPro、Za、Ta計(jì)算值與理論值之差值均小于21.1nT，第3組模型試驗(yàn)ΔTPro的計(jì)算值與理論值之差均小于4.2nT，此誤差相比于理論值可以忽略不計(jì)。這說(shuō)明只要給定一個(gè)合理的閾值δ，迭代均能收斂且能夠達(dá)到較高的計(jì)算精度。另外，轉(zhuǎn)換計(jì)算ΔTPro的誤差均低于其他分量的轉(zhuǎn)換誤差，原因在于ΔTPro與ΔT的分布形態(tài)和均值一致性更高，因而波數(shù)域奇點(diǎn)的影響也更弱。

圖7 第2組模型試驗(yàn)擬合誤差ε隨迭代次數(shù)的變化曲線

圖8 第3組模型試驗(yàn)擬合誤差ε隨迭代次數(shù)的變化曲線

圖9 第2組模型試驗(yàn)的模型理論值(上)及其與計(jì)算值之差值(下) (a)ΔTPro； (b)Za； (c)Ta； 上圖和下圖的等值線間隔分別為2000nT和0.5nT

圖10 第3組模型試驗(yàn)ΔTPro的理論值(上)及其與計(jì)算值之差(下) 主磁場(chǎng)的強(qiáng)度(T0)、傾角(I)與偏角(D)分別為：50000 nT、45°與45° (a)； 50000 nT、 90°與0° (b)； 20000 nT、90°與0° (c)； 上圖和下圖的等值線間隔分別為1000nT與0.5nT

2.2.3 計(jì)算效率

第4組模型的差異僅在于點(diǎn)距與線距的變化，旨在測(cè)試計(jì)算效率及其隨數(shù)據(jù)總量的變化特征。

對(duì)于Zhen等[12]和甄慧翔等[13-14]提出的計(jì)算方法，由于目標(biāo)函數(shù)的導(dǎo)數(shù)方程難以得到，因而采用差分代替微分計(jì)算每個(gè)測(cè)點(diǎn)處的導(dǎo)數(shù)，故計(jì)算量較大，不適用于較大數(shù)據(jù)量的情況。例如，對(duì)于101×101網(wǎng)格數(shù)據(jù)的計(jì)算，每次計(jì)算大約需要30min，對(duì)于201×201網(wǎng)格數(shù)據(jù)則需要8h。研究顯示，計(jì)算耗時(shí)隨著數(shù)據(jù)量增長(zhǎng)呈二次冪增長(zhǎng)[13]，即當(dāng)x與y方向的數(shù)據(jù)量均增大1倍時(shí)，每次迭代的耗時(shí)將增加至原來(lái)的16倍。

對(duì)于本文所提方法，第4組模型3套數(shù)據(jù)的數(shù)據(jù)量分別為801×801、401×401和201×201，在筆記本電腦(Inter(R) Xeon(R) CPU E3-1505M v6 @3.00 GHz的處理器、內(nèi)存8GB、64位操作系統(tǒng))上每次迭代分別耗時(shí)約4.0、0.9、0.2s，若設(shè)閾值δ=1nT，則3個(gè)模型的計(jì)算耗時(shí)分別為21.2、4.1、0.9s。這說(shuō)明本文算法具有較高的計(jì)算效率。

由前文理論分析可知，本文算法的計(jì)算效率主要取決于快速傅里葉正、反變換，每次迭代需要進(jìn)行1次傅里葉正變換和3次傅里葉反變換，而二維快速傅里葉變換的時(shí)間復(fù)雜度為O[M×Nlg(M×N)][31](M和N分別表示x與y方向的數(shù)據(jù)量)，因此本文所提算法的時(shí)間復(fù)雜度為O[4×imax×M×Nlg(M×N)]，其中imax為迭代總次數(shù)。

3 實(shí)際應(yīng)用

上節(jié)展示了本文算法的可靠性、穩(wěn)定性及計(jì)算效率，并分析了計(jì)算精度和迭代收斂速度隨磁場(chǎng)矢量異常的模量大小、主磁場(chǎng)大小和方向的變化關(guān)系，以及計(jì)算時(shí)間隨測(cè)點(diǎn)數(shù)量的變化關(guān)系。Zhen等[12]、甄慧翔等[13-14]、Sun等[15]和Yang等[17]已經(jīng)驗(yàn)證了在實(shí)際情況中考慮ΔT近似為ΔTPro誤差的必要性。因此，此節(jié)僅討論本文算法在實(shí)際應(yīng)用中需要注意的問(wèn)題，及其與傳統(tǒng)方法計(jì)算結(jié)果之間的差異。

研究區(qū)位于河北省東部的麻城鐵礦區(qū)，該鐵礦為大型隱伏條帶狀鐵建造(Banded Iron Formation, BIF)礦體[32]。圖11a為研究區(qū)部分地面實(shí)測(cè)總磁場(chǎng)強(qiáng)度異常分布，測(cè)點(diǎn)點(diǎn)距和線距均為25m，數(shù)據(jù)點(diǎn)數(shù)為181(北向)×129(東向)。由圖可見(jiàn)，總磁異常幅值范圍為-2806～11402nT。當(dāng)?shù)氐闹鞔艌?chǎng)強(qiáng)度為54000nT，傾角和偏角分別為57.9°和-7.1°[16]。根據(jù)文獻(xiàn)[8-9]，采用投影分量ΔT近似ΔTPro的最大誤差Emax=Ta2/(2T0)，可以估算出研究區(qū)該最大誤差約為1204nT。如此大的計(jì)算誤差將對(duì)后續(xù)的定量解釋產(chǎn)生比較嚴(yán)重的影響，因此必須考慮投影分量ΔTPro與實(shí)測(cè)ΔT之間的差異。

為了壓制快速傅里葉變換的邊界效應(yīng)，并解決位場(chǎng)數(shù)據(jù)波數(shù)域轉(zhuǎn)換處理的奇點(diǎn)問(wèn)題，需要對(duì)實(shí)測(cè)ΔT進(jìn)行背景場(chǎng)去除及擴(kuò)邊處理。從圖11a可以看出，研究區(qū)域主要具有兩個(gè)明顯的磁力異常：中北部的正值異常范圍較小，且在其北側(cè)具有明顯的負(fù)值異常旁瓣；中南部的正值異常幅值較高、范圍較大，沒(méi)有明顯的負(fù)值異常旁瓣。由于這兩個(gè)磁異常相距較近，推測(cè)其場(chǎng)源磁化方向相同，因此猜測(cè)測(cè)區(qū)北部大面積的負(fù)值異?？赡懿糠质悄喜看竺娣e正值異常的旁瓣或延伸。由于該測(cè)區(qū)范圍較小、背景場(chǎng)并不明顯，直接地去除背景場(chǎng)極有可能使得目標(biāo)磁力異常發(fā)生畸變，因此不進(jìn)行背景場(chǎng)去除而僅進(jìn)行擴(kuò)邊處理。根據(jù)測(cè)區(qū)南部大面積正值異常的空間分布特征，沿南北向和東西向均進(jìn)行擴(kuò)邊，選擇兩種擴(kuò)邊方案：500m和1000m。擴(kuò)邊后的總數(shù)據(jù)量分別為221×169、261×209，擴(kuò)邊處理之后的結(jié)果如圖11b和圖11c所示，參數(shù)統(tǒng)計(jì)見(jiàn)表3。另外，在研究區(qū)的西北邊界還存在一個(gè)幅值略高的正值異常，該異常并不完整，其轉(zhuǎn)換計(jì)算誤差較大，所以本文不予關(guān)注。

圖11 研究區(qū)域?qū)崪y(cè)總磁場(chǎng)強(qiáng)度異常(ΔT)及擴(kuò)邊結(jié)果 (a)實(shí)測(cè)數(shù)據(jù)； (b)四周均擴(kuò)邊500m； (c)四周均擴(kuò)邊1000m 等值線間距為1000 nT； 紅色、藍(lán)色與黑色邊框分別表示未擴(kuò)邊及四周均擴(kuò)邊500m、1000m后的區(qū)域邊界

采用本文算法，對(duì)圖11所示3套數(shù)據(jù)分別進(jìn)行計(jì)算。擬合誤差ε隨迭代次數(shù)的變化曲線見(jiàn)圖12，可見(jiàn)擴(kuò)邊范圍越大，最終迭代收斂時(shí)的擬合差越小，但擬合差降低趨勢(shì)越來(lái)越不明顯。說(shuō)明在實(shí)際數(shù)據(jù)處理時(shí)，適當(dāng)?shù)臄U(kuò)邊是必要的，過(guò)度擴(kuò)邊會(huì)極大地增加計(jì)算量。從圖12還看出，本文方法迭代收斂速度較快、計(jì)算穩(wěn)定且計(jì)算精度較高。

圖13為采用本文方法、5次迭代后的ΔT擬合差分布，其參數(shù)統(tǒng)計(jì)見(jiàn)表4?？梢?jiàn)，在未擴(kuò)邊時(shí)，ΔT擬合差幅值較大；而擴(kuò)邊之后，ΔT擬合差的平均值快速降低至幾個(gè)nT量級(jí)；擴(kuò)邊與不擴(kuò)邊情況下ΔT擬合差的幅值變化范圍與標(biāo)準(zhǔn)差幾乎一致，說(shuō)明通過(guò)擴(kuò)邊可以有效改善磁力異常波數(shù)域轉(zhuǎn)換計(jì)算的奇點(diǎn)問(wèn)題。

圖14是基于擴(kuò)邊1000m、迭代5次計(jì)算得到的Ha,x、Ha,y、Za、Ta分布。由圖可知，前文關(guān)于研究區(qū)存在兩個(gè)較大的ΔT磁力異常及其場(chǎng)源的推斷是合理的。圖15是采用本文方法與直接采用傳統(tǒng)波數(shù)域方法計(jì)算結(jié)果(ΔTPro、Ha,x、Ha,y、Za、Ta)的差值，其數(shù)據(jù)統(tǒng)計(jì)結(jié)果見(jiàn)表5?？梢?jiàn)兩種方法的計(jì)算結(jié)果差值最高可達(dá)600nT，且該差異與磁異常在空間分布上關(guān)系密切。若忽略ΔT近似為ΔTPro的誤差，直接進(jìn)行分量與模量轉(zhuǎn)換，會(huì)產(chǎn)生較大的誤差，該誤差對(duì)于磁異常的定性解釋影響較小，但對(duì)于定量解釋的影響較大，這再次說(shuō)明在高磁環(huán)境中考慮ΔT近似為投影分量ΔTPro的誤差是非常必要的。此外，由圖15a可以發(fā)現(xiàn)，ΔTPro均小于ΔT，這與前文理論分析結(jié)論一致，即ΔT≥ΔTPro，這也從另一個(gè)側(cè)面證明了本文方法的可靠性。

圖12 在不同擴(kuò)邊大小情況下擬合誤差ε隨迭代次數(shù)的變化曲線

表3 實(shí)測(cè)ΔT磁異常數(shù)據(jù)擴(kuò)邊結(jié)果統(tǒng)計(jì)表

表4 本文方法計(jì)算的ΔT擬合差統(tǒng)計(jì)參數(shù)表

圖13 研究區(qū)采用本文方法計(jì)算的ΔT擬合差 (a)未擴(kuò)邊； (b)四周均擴(kuò)邊500m； (c)四周均擴(kuò)邊1000m

圖14 研究區(qū)數(shù)據(jù)采用本文方法的計(jì)算結(jié)果 (a)ΔTPro； (b)Ha,x； (c)Ha,y； (d)Za； (e)Ta 等值線間距為1000nT

表5 本文方法與傳統(tǒng)方法計(jì)算結(jié)果及其差值統(tǒng)計(jì)

圖15 研究區(qū)采用本文方法與傳統(tǒng)波數(shù)域轉(zhuǎn)換方法計(jì)算結(jié)果的差值 (a)ΔTPro； (b)Ha,x； (c)Ha,y； (d)Za； (e)Ta 等值線間距為50nT

4 結(jié)論

在高磁環(huán)境中，將ΔT近似為ΔTPro的誤差較大，若忽略該誤差，對(duì)實(shí)測(cè)ΔT的數(shù)據(jù)預(yù)處理、正演、反演及定量解釋會(huì)產(chǎn)生較大的影響。雖然已有學(xué)者提出了基于ΔT精確計(jì)算ΔTPro的方法，但是計(jì)算量大、計(jì)算效率低，限制了實(shí)際應(yīng)用。為此，本文提出了一種高效率的迭代計(jì)算方法，通過(guò)理論分析、模型試驗(yàn)與實(shí)際應(yīng)用，得到如下結(jié)論。

(1)本文方法具有較快的迭代收斂速度、較高的計(jì)算穩(wěn)定性、計(jì)算精度及計(jì)算效率，一般僅需2～3次迭代即可將ΔT擬合差的幅值控制在幾個(gè)nT之內(nèi)。

(2)在每次迭代計(jì)算ΔT近似為ΔTPro的誤差時(shí)，采用式(8)或采用實(shí)測(cè)ΔT數(shù)據(jù)計(jì)算ΔT(i)，迭代的收斂性沒(méi)有顯著差異。但是本文作者認(rèn)為式(5)所示誤差E與ΔT、T0、Ta之間是耦合的，應(yīng)采用式(8)計(jì)算；若采用實(shí)測(cè)ΔT數(shù)據(jù)代替可能會(huì)致使迭代過(guò)程存在振蕩的風(fēng)險(xiǎn)。

(3)采用相鄰兩次迭代的擬合誤差之差δ作為迭代終止條件比采用擬合誤差更合理。

(4)本文方法的迭代收斂速度主要受控于磁場(chǎng)矢量異常的模量大小Ta，而受主磁場(chǎng)強(qiáng)度的影響較弱。

(5)本文方法的計(jì)算效率主要取決于快速傅里葉正、反變換，每次迭代需要進(jìn)行1次正變換和3次反變換，因此本文所提算法的時(shí)間復(fù)雜度為O[4×imax×M×Nlg(M×N)]。

(6)本文方法的計(jì)算誤差主要來(lái)源于快速傅里葉正、反變換的計(jì)算精度，其根源在于邊界效應(yīng)與磁力異常波數(shù)域轉(zhuǎn)換計(jì)算的奇點(diǎn)問(wèn)題，背景場(chǎng)、模量Ta與主磁場(chǎng)方向或磁化方向?qū)τ?jì)算精度的影響均表現(xiàn)于此，對(duì)研究區(qū)進(jìn)行適當(dāng)擴(kuò)邊可以有效降低計(jì)算誤差。

(7)背景場(chǎng)對(duì)波數(shù)域轉(zhuǎn)換存在較大影響，但是“粗暴”地去除背景場(chǎng)極有可能導(dǎo)致目標(biāo)磁力異常發(fā)生畸變，因此在實(shí)際應(yīng)用時(shí)應(yīng)該事先仔細(xì)分析ΔT的分布特征，根據(jù)研究目標(biāo)選擇是否剔除背景場(chǎng)。另外，適當(dāng)?shù)臄U(kuò)邊本身也可以壓制背景場(chǎng)效應(yīng)。

(8)若忽略ΔT近似為ΔTPro的誤差而直接進(jìn)行分量與模量轉(zhuǎn)換，會(huì)產(chǎn)生較大的誤差，該差異對(duì)于磁異常定性解釋影響較小，但對(duì)于定量解釋的影響較大。因此，在高磁環(huán)境中，考慮ΔT近似為投影分量ΔTPro的誤差是非常必要的。

筆者已經(jīng)將本文算法編寫為計(jì)算機(jī)軟件，讀者可與本文作者聯(lián)系免費(fèi)獲取該軟件程序。

本文所用實(shí)際地面磁測(cè)數(shù)據(jù)來(lái)源于中國(guó)冶金地質(zhì)總局礦產(chǎn)資源研究院，在此表示衷心的感謝！