基于磁異常ΔT計算投影分量ΔTPro的迭代算法

胡正旺 杜勁松*③ 孫石達 陳 超

(①中國地質大學(武漢)地球物理與空間信息學院，湖北武漢 430074；②地球內部多尺度成像湖北省重點實驗室，湖北武漢 430074；③地質過程與礦產資源國家重點實驗室，湖北武漢 430074；④河北地質大學勘查技術與工程學院，河北石家莊 050031)

0 引言

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

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

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

1 方法原理與算法流程

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

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

T=T0+Ta

(1)

如引言所述，實測總磁場強度異常ΔT往往表示為總磁場強度|T|與主磁場強度|T0|(即T0)之差

ΔT=|T|-|T0|

(2)

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

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

ΔTPro=|Ta|cosθ

(3)

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

E=ΔT-ΔTPro

(4)

根據(jù)圖1，經過推導與化簡可得誤差E的表達式[8]為

(5)

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

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

(6)

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

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

(7)

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

(8)

(9)

(10)

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

(11)

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

(12)

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

(1)對實測數(shù)據(jù)ΔT進行數(shù)據(jù)預處理，依次進行奇異值剔除、數(shù)據(jù)網格化、擴邊、初步去噪；

(2)給定ε0；

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

2 模型試驗

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

2.1 理論模型

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

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

表1 理論磁性模型和磁場觀測參數(shù)表

2.2 試驗結果及其分析

2.2.1 方法有效性

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2.3 計算效率

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

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

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

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

3 實際應用

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

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

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

圖11 研究區(qū)域實測總磁場強度異常(ΔT)及擴邊結果 (a)實測數(shù)據(jù)； (b)四周均擴邊500m； (c)四周均擴邊1000m 等值線間距為1000 nT； 紅色、藍色與黑色邊框分別表示未擴邊及四周均擴邊500m、1000m后的區(qū)域邊界

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

圖13為采用本文方法、5次迭代后的ΔT擬合差分布，其參數(shù)統(tǒng)計見表4?？梢姡谖磾U邊時，ΔT擬合差幅值較大；而擴邊之后，ΔT擬合差的平均值快速降低至幾個nT量級；擴邊與不擴邊情況下ΔT擬合差的幅值變化范圍與標準差幾乎一致，說明通過擴邊可以有效改善磁力異常波數(shù)域轉換計算的奇點問題。

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

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

表3 實測ΔT磁異常數(shù)據(jù)擴邊結果統(tǒng)計表

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

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

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

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

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

4 結論

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

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

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

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

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

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

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

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

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

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

本文所用實際地面磁測數(shù)據(jù)來源于中國冶金地質總局礦產資源研究院，在此表示衷心的感謝！