改進(jìn)的低緯度化極穩(wěn)定算法

張 琪 張英堂 李志寧 范紅波(陸軍工程大學(xué)石家莊校區(qū)車(chē)輛與電氣工程系,河北石家莊 050003； ②94019部隊(duì),新疆和田 848000)

1 引言

化極，即把觀測(cè)到的磁力異常歸算為假定磁源體位于地磁極處產(chǎn)生的磁力異常。頻率域常規(guī)化極因子(Routine RTP Operator,RO)屬于放大類(lèi)因子，在低緯度地區(qū)，RO的幅值變得極大，導(dǎo)致化極極不穩(wěn)定，會(huì)沿著地磁偏角方向產(chǎn)生嚴(yán)重的條帶狀異常。在噪聲環(huán)境下，RO對(duì)高頻成分的放大作用較強(qiáng)，導(dǎo)致化極結(jié)果的信噪比較低。因此，有必要研究快速、穩(wěn)定的低緯度化極方法。

低緯度化極分為空間域方法和頻率域方法。Baranov[1]首先提出化極的概念，并將空間域化極視為一種褶積運(yùn)算；后來(lái)，空間域化極又被視為空間域反演問(wèn)題[2,3]?？臻g域方法的優(yōu)點(diǎn)是精度高，缺點(diǎn)是計(jì)算復(fù)雜且計(jì)算量大，涉及到大型線性方程組的求解問(wèn)題。Bhattacharyya[4]利用傅里葉變換將空間域復(fù)雜的褶積運(yùn)算轉(zhuǎn)換成頻率域簡(jiǎn)單的乘積運(yùn)算，實(shí)現(xiàn)了頻率域化極。頻率域化極方法原理簡(jiǎn)單、計(jì)算快捷，已成為低緯度地區(qū)磁異常定性解釋的有效手段[5]。根據(jù)壓制區(qū)間的不同，可將頻率域化極方法劃分為兩類(lèi)： 一類(lèi)主要用來(lái)壓制不穩(wěn)定的放大區(qū)，即磁偏角的垂直方向及其附近區(qū)域，稱(chēng)之為“局部壓制類(lèi)”化極方法； 另一類(lèi)主要從全局進(jìn)行壓制，稱(chēng)之為“全局壓制類(lèi)”化極方法。

局部壓制類(lèi)化極方法主要對(duì)放大區(qū)內(nèi)的RO進(jìn)行改進(jìn)，按改進(jìn)方式的不同又可分為濾波法、阻尼法和分類(lèi)討論法等三種。

濾波法指在放大區(qū)將RO乘以一個(gè)特殊的濾波因子。姚長(zhǎng)利等[6]利用一個(gè)余弦函數(shù)濾波因子定義壓制因子(Suppression Filter,SF)，但壓制邊界角較難確定； Li[7]提出方位角濾波法(Azimuthal Filtering,AF)，其實(shí)質(zhì)是給RO乘以一個(gè)正弦函數(shù)濾波因子。

阻尼法指在放大區(qū)將RO的分母加上一個(gè)較小的阻尼因子。姚長(zhǎng)利等[8]定義了直接阻尼(Direct Damping,DD)法，并利用一個(gè)余弦阻尼因子壓制不穩(wěn)定的放大區(qū)； 林曉星等[9]通過(guò)增加相位因子和頻率因子改進(jìn)DD法，并將其定義為變頻雙向阻尼因子(Frequency Conversion Bidirectional Damping Factor，F(xiàn)C)法。

分類(lèi)討論法包括反對(duì)稱(chēng)因子(Antisymmetric Factor，ASF)法和非線性閾值法(Nonlinear Thresholding Method,NT)。Guo等[10]利用一個(gè)角度閾值將化極區(qū)域分成三個(gè)不同區(qū)域，在放大區(qū)內(nèi)應(yīng)用反對(duì)稱(chēng)因子(ASF)，在其他區(qū)域運(yùn)用RO以保留有效的化極特征； Zhang等[11]把RO表示為復(fù)數(shù)形式，并用一個(gè)最優(yōu)非線性幅度閾值分別穩(wěn)定RO的實(shí)部和虛部。

全局壓制類(lèi)化極方法的基本原理是利用數(shù)學(xué)方程構(gòu)造新的化極因子以穩(wěn)定全部化極區(qū)域。Mendon?a等[12]在向上延拓的基礎(chǔ)上，推導(dǎo)了截?cái)嗟奶├占?jí)數(shù)法(Truncated Taylor Series Approximation Method,TT)，但該法的參數(shù)較難確定；雙曲正弦函數(shù)法(Hyperbolic Sine Method,HS)[13]是將RO的分子和分母分別減去和加上一個(gè)改進(jìn)的雙曲正弦函數(shù)來(lái)穩(wěn)定RO。

學(xué)者們針對(duì)低緯度化極做了深入研究，提出許多頻率域化極方法，但仍存在以下問(wèn)題： ①RO在垂直于地磁偏角方向極不穩(wěn)定甚至奇異，導(dǎo)致化極結(jié)果出現(xiàn)嚴(yán)重的條帶狀干擾； ②局部壓制類(lèi)和全局壓制類(lèi)化極方法通過(guò)改進(jìn)RO壓制放大區(qū)，該過(guò)程不可避免地引入?yún)?shù)，并且參數(shù)的取值直接影響化極結(jié)果的精度，但針對(duì)參數(shù)取值的研究較少，往往都是利用經(jīng)驗(yàn)值或給定參數(shù)范圍； ③參數(shù)最優(yōu)的RO在較低緯度仍表現(xiàn)出明顯的高通濾波特性，對(duì)高頻成分具有較強(qiáng)的放大作用，并且緯度越低，作用越明顯。

針對(duì)以上問(wèn)題，對(duì)RO進(jìn)行改造，提出了局部壓制類(lèi)和全局壓制類(lèi)化極方法，然后通過(guò)計(jì)算化極結(jié)果與歸一化磁源強(qiáng)度的最大相關(guān)系數(shù)得到了最優(yōu)參數(shù)和最優(yōu)參數(shù)下的化極結(jié)果，最后提出了正則化改進(jìn)算法，得到了穩(wěn)定的低緯度化極結(jié)果。

2 化極理論

2.1 頻率域化極原理

頻率域總場(chǎng)強(qiáng)度的RO可表示為

(1)

(2)

式中θ=arctan(v/u)。H(θ)依賴(lài)于I和D。在低緯度地區(qū)，I的值趨近于零，RO將變得不穩(wěn)定。在極端情況下，即當(dāng)θ=D±90°時(shí)，H(θ)的值趨近于無(wú)窮。設(shè)β為一小角度，則H(θ)在θ±β的區(qū)域內(nèi)放大作用最強(qiáng)，稱(chēng)之為“放大區(qū)”。

2.2 局部壓制類(lèi)方法

濾波法的基本原理是在放大區(qū)將RO乘以一個(gè)濾波因子，在其他區(qū)域仍然使用RO，其原理表示為

(3)

式中：Hf(θ)表示濾波法的化極因子；Fβ(θ)表示濾波因子；α=θ-θ0,θ0=D±90°。

定義SF的余弦濾波因子和AF的正弦濾波因子分別為FSF(θ)和FAF(θ)，其表達(dá)式為

(4)

式中：Iac為低緯度特征角；p為壓制衰減率。

阻尼法的實(shí)質(zhì)是在放大區(qū)將RO的分母加上一個(gè)阻尼因子，其原理可表示為

Hd(θ)

(5)

式中：Hd(θ)表示阻尼法的化極因子；Fd(θ)表示阻尼因子。

用FDD(θ)和FFC(θ)分別表示DD法和FC法的阻尼因子，其表達(dá)式為

(6)

式中d、k和h分別表示阻尼因子、相位因子和頻率因子。

用HASF(θ)表示ASF法的化極因子，其表達(dá)式為

(7)

式中θ′表示角度閾值。NT的基本原理可描述為

(8)

式中：x表示RO的虛部或?qū)嵅浚粁′表示新的虛部或?qū)嵅?；?和θ2分別表示放大區(qū)的兩個(gè)邊緣角；Q表示非線性幅值閾值。

2.3 全局壓制類(lèi)方法

設(shè)HTT(θ)表示TT的化極因子，其定義為

HTT(θ)= {[sinI-icosIcos(θ-D)]×

(9)

式中：M為一正整數(shù)，表示泰勒級(jí)數(shù)展開(kāi)項(xiàng)的系數(shù)；K=2π|r|；l表示相對(duì)于觀測(cè)面的延拓距離。

用HHS(θ)表示HS的化極因子，其表達(dá)式為

HHS(θ)=

(10)

式中：sinh表示雙曲正弦函數(shù)；參數(shù)A的取值范圍為2≤A≤12；參數(shù)B的取值范圍為0.0001≤B≤0.1。

3 尋優(yōu)

3.1 化極結(jié)果尋優(yōu)

圖1所示為頻率域化極算法流程。根據(jù)圖1可知，磁異?；瘶O過(guò)程中，首先獲得待測(cè)區(qū)域的磁異常信息，然后分別應(yīng)用常規(guī)化極方法、局部壓制類(lèi)化極方法和全局壓制類(lèi)化極方法對(duì)總場(chǎng)數(shù)據(jù)進(jìn)行化極，并利用相關(guān)系數(shù)法對(duì)化極因子的參數(shù)尋優(yōu)，計(jì)算最優(yōu)參數(shù)下的化極結(jié)果，并對(duì)化極效果進(jìn)行評(píng)判。若化極結(jié)果能夠準(zhǔn)確、清晰地突出磁源體的磁異常特征，表明化極效果較理想，即可直接輸出最優(yōu)化極結(jié)果；若化極效果不理想，尤其當(dāng)緯度很低且噪聲很大時(shí)，參數(shù)最優(yōu)的化極因子仍存在明顯的高頻放大效應(yīng)，此時(shí)需對(duì)化極因子進(jìn)行正則化改進(jìn)，計(jì)算并輸出正則化改進(jìn)后的化極結(jié)果。

圖1 化極算法流程

3.2 相關(guān)系數(shù)法

歸一化磁源強(qiáng)度(Normalized Source Strength,NSS)是由磁偶極子的磁梯度張量的特征值矩陣推導(dǎo)出的一個(gè)旋轉(zhuǎn)不變量[14,15]，NSS弱敏感于磁化方向，即使在低緯度地區(qū)，NSS仍能反映磁源體的磁性特征。

局部壓制類(lèi)和全局壓制類(lèi)化極方法可對(duì)不穩(wěn)定的放大區(qū)進(jìn)行有效壓制，卻不可避免地引入了參數(shù)，而參數(shù)取值直接影響化極結(jié)果的精度。因此，為了提高化極結(jié)果的精度，利用相關(guān)系數(shù)法對(duì)化極因子的參數(shù)尋優(yōu)。相關(guān)系數(shù)可以描述化極結(jié)果與NSS的相關(guān)程度。相關(guān)系數(shù)越大，表明二者的對(duì)應(yīng)關(guān)系越好，即最大相關(guān)系數(shù)對(duì)應(yīng)最優(yōu)參數(shù)。相關(guān)系數(shù)計(jì)算公式為

(11)

式中：X表示化極結(jié)果； cov(·，·)為協(xié)方差算子；D(·)為方差算子。

3.3 Tikhonov正則化

頻率域化極可表示為

ST2(u,v)=H(u,v)ST1(u,v)

(12)

式中：ST1(u,v)表示斜磁化的總場(chǎng)數(shù)據(jù)的頻譜；ST2(x,y)表示化極結(jié)果的頻譜。

在低緯度地區(qū)，H(u,v)表現(xiàn)出明顯的高通濾波特性。Zeng等[16]和Yin等[17]采用正則化濾波算子消除磁異常轉(zhuǎn)換因子的噪聲放大效應(yīng)，有效提高了磁異常的信噪比。由于式(12)為病態(tài)問(wèn)題，可將其轉(zhuǎn)化為最小二乘形式[18,19]

min{‖H-1(u,v)ST2(u,v)-ST1(u,v)‖22+

λ‖ST2(u,v)‖22}

(13)

式中： ‖·‖22表示歐幾里德范數(shù)的平方；λ為正則化參數(shù)。

求解式(13)可得

(14)

若Hλ(u,v)表示化極因子的正則化改進(jìn)形式，可表示為

(15)

式(15)的極坐標(biāo)形式可表示為

(16)

式中H′(θ)表示任意參數(shù)最優(yōu)的化極因子。

采用L—曲線法[20]對(duì)λ進(jìn)行參數(shù)尋優(yōu)。L—曲線是一個(gè)方便的雙對(duì)數(shù)圖，在L—曲線的拐點(diǎn)處，即曲率最大的位置，對(duì)應(yīng)最優(yōu)λ，此時(shí)化極因子的壓制效果最好。

3.4 均方根誤差

為了衡量化極結(jié)果的精度，引入均方根誤差(RMSE)作為定量評(píng)價(jià)指標(biāo)，其定義為

(17)

式中：Nx和Ny分別表示x和y方向上的采樣點(diǎn)總數(shù)；X⊥(i,j)表示垂直磁化條件下的理論磁異常值；X(i,j)表示化極磁異常值。

根據(jù)RMSE的定義，理論上，RMSE的值越小，則化極結(jié)果越接近于垂直磁化條件下的磁異常值，化極結(jié)果的精度越高；RMSE的值越大，表明化極結(jié)果的精度越低。

4 仿真試驗(yàn)

假設(shè)空間中存在如圖2所示的由長(zhǎng)方體、正方體和圓柱體組成的磁源體模型。長(zhǎng)方體、正方體和圓柱體的中心坐標(biāo)分別為(30m,-10m,16m)、(10m,30m,14m)和(-25m,-5m,18m)，長(zhǎng)方體和正方體的邊長(zhǎng)分別為(20m,40m,20m)和(20m,20m,20m)，圓柱體的半徑和母線長(zhǎng)度分別為10m和50m。模型網(wǎng)格數(shù)為121×121，水平方向的采樣間隔為1m。

圖2 組合體的空間位置示意圖

由于地磁偏角僅起到圖形旋轉(zhuǎn)的作用，因此假設(shè)D=23°，組合體的磁化強(qiáng)度大小均為30A/m。對(duì)總場(chǎng)強(qiáng)度數(shù)據(jù)添加均值為0、方差為100nT的高斯白噪聲，并且所有仿真均在該噪聲環(huán)境下進(jìn)行。需要說(shuō)明的是，仿真中的地磁傾角、地磁偏角、磁化強(qiáng)度和噪聲大小等參數(shù)均為任意取值，也可以選擇其他數(shù)值驗(yàn)證所述方法的正確性。

當(dāng)I分別等于90°、10°和1°時(shí)，組合體在z=0平面的總場(chǎng)強(qiáng)度如圖3所示。根據(jù)圖3可知，垂直磁化條件下的磁異常(圖3a)與組合體的水平位置的對(duì)應(yīng)關(guān)系最佳，磁異常形態(tài)清晰，易于辨識(shí)。相比之下，受磁化方向的影響，低緯度地區(qū)的磁異常(圖3b和圖3c)形態(tài)較復(fù)雜，與組合體的對(duì)應(yīng)關(guān)系較差，磁異常數(shù)據(jù)的解釋較困難。因此，為了降低斜磁化的影響，有必要對(duì)低緯度磁異常數(shù)據(jù)進(jìn)行化極。

當(dāng)I=10°時(shí)，利用前述九種方法進(jìn)行化極的結(jié)果如圖4所示，其RMSE如表1所示。

由圖4可知，常規(guī)化極結(jié)果(圖4i)產(chǎn)生了較多的沿著地磁偏角方向的條帶狀異常，化極效果不理想。相比之下，局部壓制類(lèi)和全局壓制類(lèi)化極方法的化極結(jié)果(圖4a～圖4h)的條帶狀異常明顯減少，磁異常形態(tài)較清晰，辨識(shí)度較高。由表1可知，與RO法相比，SF、DD、ASF、AF、FC、NT、TT和HS法結(jié)果的RMSE較小，表明局部壓制類(lèi)和全局壓制類(lèi)化極方法對(duì)放大區(qū)的壓制效果較好，化極結(jié)果的精度較高。根據(jù)上述分析可知，當(dāng)I=10°時(shí)，局部壓制類(lèi)和全局壓制類(lèi)化極方法的化極效果優(yōu)于常規(guī)化極方法，基本達(dá)到低緯度化極的目的。

當(dāng)I=1°時(shí)，九種不同方法的化極結(jié)果如圖5所示，其RMSE如表2所示。

由圖5可知，ASF(圖5c)和常規(guī)化極(圖5i)結(jié)果已被噪聲完全掩蓋，磁異常值被嚴(yán)重放大，磁異常形態(tài)難以辨識(shí)；SF(圖5a)、DD(圖5b)和FC(圖5e)化極結(jié)果產(chǎn)生了大量沿地磁偏角方向的條帶狀干擾，辨識(shí)度較差；AF(圖5d)和NT(圖5f)的化極結(jié)果產(chǎn)生了較多的條帶狀異常；相比之下，TT(圖5g)和HS(圖5h)化極結(jié)果的條帶狀異常基本沒(méi)有出現(xiàn)，磁異常形態(tài)較清晰。由表2可知，RO和ASF計(jì)算結(jié)果的RMSE非常大，表明這兩種算法的化極結(jié)果精度最低、化極效果差；與RO和ASF相比，SF、DD、AF、FC、NT、TT和HS的RMSE較小，表明這七種方法的化極精度較高。根據(jù)上述分析，當(dāng)I=1°時(shí)，全局壓制類(lèi)化極方法的整體化極效果要明顯優(yōu)于局部壓制類(lèi)方法和常規(guī)化極方法，全局壓制類(lèi)化極方法基本能夠達(dá)到低緯度化極的目的。

圖3 組合體在不同地磁傾角時(shí)的總磁場(chǎng)強(qiáng)度分布圖(z=0)(a)I=90°； (b)I=10°； (c)I=1°

圖4 九種方法在I=10°時(shí)的磁異?；瘶O結(jié)果

(a)SF：β=20°，Corr=0.761； (b)DD：d=0.09，Corr=0.762； (c)ASF：θ′=64°，Corr=0.763； (d)AF：p=1，Corr=0.759； (e)FC：k=1，h=2，Corr=0.763； (f)NT：Q=15，Corr=0.759； (g)TT：M=20，Corr=0.755； (h)HS：A=5，B=0.001，Corr=0.718； (i)RO：Corr=0.756

表1 I=10°時(shí)九種方法磁異常化極結(jié)果的RMSE

圖5 九種方法在I=1°時(shí)的化極結(jié)果(a)SF：β=20°，Corr=0.649； (b)DD：d=0.5，Corr=0.72； (c)ASF：θ′=80°，Corr=0.216； (d)AF：p=5，Corr=0.736； (e)FC：k=2，h=2，Corr=0.735； (f)NT：Q=3，Corr=0.749； (g)TT：M1=20，Corr=0.755； (h)HS：A=12，B=0.03，Corr=0.745； (i)RO：Corr=0.131

表2 I=1°時(shí)九種方法的磁異?；瘶O結(jié)果的RMSE

為了分析圖5中九種化極結(jié)果存在明顯差異的原因，分別繪制九種化極因子在I=1°時(shí)的幅頻特性曲線(圖6)。

由圖6可知，當(dāng)I=1°時(shí)，在不穩(wěn)定區(qū)域化極因子幅值水平RO(圖6i)最高，ASF(圖6c)次之，然后依次是SF(圖6a)、DD(圖6b)、FC(圖6e)、AF(圖6d)、NT(圖6f)，TT(圖6g)和HS(圖6h)最低。綜合圖5和圖6可知，化極因子在不穩(wěn)定區(qū)域內(nèi)的幅值水平越低，化極因子對(duì)高頻成分的壓制作用越強(qiáng)，化極就越穩(wěn)定。

當(dāng)I=1°時(shí)，為了有效抑制化極因子的高頻放大效應(yīng)，利用式(16)對(duì)其進(jìn)行正則化改進(jìn)。正則化改進(jìn)后的九種化極結(jié)果如圖7所示，其RMSE統(tǒng)計(jì)結(jié)果如表3所示、幅頻特性曲線如圖8所示。

對(duì)比圖7與圖5可以看出，SF(圖7a)、DD(圖7b)、ASF(圖7c)、AF(圖7d)、FC(圖7e)和常規(guī)化極結(jié)果(圖7i)均得到明顯改善，條帶狀異?；鞠?，磁異常形態(tài)清晰、辨識(shí)度較高；NT(圖7f)和HS(圖7h)化極結(jié)果的分辨率更高；而TT(圖7g)的化極結(jié)果無(wú)明顯變化。

圖6 I=1°時(shí)九種化極因子的幅頻特性曲線(a)SF； (b)DD； (c)ASF； (d)AF； (e)FC； (f)NT； (g)TT； (h)HS； (i)RO

圖7 正則化改進(jìn)后九種方法在I=1°時(shí)的化極結(jié)果(a)SF； (b)DD； (c)AS； (d)AF； (e)FC； (f)NT； (g)TT； (h)HS； (i)RO

表3與表2相比，SF、DD、ASF、AF、FC、HS和RO的RMSE均明顯變小，其中，RO和ASF的RMSE下降幅度較大；NT的RMSE增大，而TT的RMSE基本不變。上述分析表明，正則化改進(jìn)后，SF、DD、AF、ASF、FC、HS和RO的化極精度得到明顯提高。需要注意的是，NT的RMSE增大是由于正則化改進(jìn)后的化極結(jié)果產(chǎn)生了一定的偏置，導(dǎo)致磁異常數(shù)值整體變小。由圖7可知，NT(圖7f)的化極結(jié)果能夠清晰地反映出組合體的磁異常分布，化極效果較理想，基本達(dá)到化極的目的。

表3 I=1°時(shí)正則化改進(jìn)后九種方法的磁異?；瘶O結(jié)果的RMSE

圖8 I=1°時(shí)正則化改進(jìn)后九種化極因子的幅頻特性曲線(a)SF； (b)DD； (c)ASF； (d)AF； (e)FC； (f)NT； (g)TT； (h)HS； (i)RO

對(duì)比圖8與圖6可以看出，SF(圖8a)、DD(圖8b)、ASF(圖8c)、AF(圖8d)、FC(圖8e)、NT(圖8f)和RO(圖8i)在不穩(wěn)定區(qū)域內(nèi)的幅值均明顯降低，其中，ASF和RO的降低幅度較大；TT和HS無(wú)明顯變化。上述幅值變化對(duì)應(yīng)在圖7中表現(xiàn)為化極結(jié)果的條帶狀異常得到消除，磁異常形態(tài)更加清晰。

根據(jù)上述分析可知，正則化改進(jìn)算法通過(guò)降低化極因子在不穩(wěn)定區(qū)域的幅值實(shí)現(xiàn)對(duì)高頻成分的有效壓制，表現(xiàn)為化極結(jié)果的條帶狀異常基本被消除，磁異常形態(tài)更清晰，化極更加穩(wěn)定，化極結(jié)果的精度得到提高。

5 實(shí)驗(yàn)

為了驗(yàn)證改進(jìn)化極算法的實(shí)際應(yīng)用效果，對(duì)巴西中部的Goiás州北部某礦區(qū)的磁異常數(shù)據(jù)[11]進(jìn)行處理。網(wǎng)格數(shù)為181×181，水平方向的采樣間隔為250m。礦區(qū)的地磁傾角和地磁偏角分別為-5°和-15°，其實(shí)測(cè)總場(chǎng)強(qiáng)度如圖9所示。

圖9 巴西中部某礦區(qū)的實(shí)測(cè)磁異常分布圖

由圖9可知，實(shí)測(cè)磁異常不僅存在嚴(yán)重的噪聲干擾，而且由于受到低緯度的影響，其磁異常形態(tài)十分復(fù)雜，磁異常值大多為負(fù)值，給該區(qū)域的磁異常數(shù)據(jù)解釋帶來(lái)了較大困難。

圖10為最優(yōu)參數(shù)下的九種方法化極結(jié)果。由圖10可知，SF(圖10a)、DD(圖10b)、ASF(圖10c)、AF(圖10d)、FC(圖10e)和常規(guī)化極結(jié)果(圖10i)均被高頻噪聲完全掩蓋，磁異常值被嚴(yán)重放大，并且產(chǎn)生了大量的條帶狀干擾，辨識(shí)度較差；NT(圖10f)的化極結(jié)果產(chǎn)生了較多的條帶狀干擾，但尚能識(shí)別部分磁異常特征； 相比之下， TT(圖10g)和HS(圖10h)中的條帶狀異常較少，礦區(qū)的正負(fù)磁異常差異較明顯，化極效果較為理想。

為了消除SF、DD、ASF、AF、FC、NT和RO這七種化極因子存在的高頻放大效應(yīng)，對(duì)其進(jìn)行正則化改進(jìn)，而不對(duì)TT和HS作任何處理。改進(jìn)后的化極結(jié)果如圖11所示，成圖時(shí)采用相同的色標(biāo)范圍。其中TT(圖11g)和HS(圖11h)為最優(yōu)參數(shù)下的化極結(jié)果，僅用于對(duì)比。對(duì)比圖11與圖10可以看出，SF(圖11a)、DD(圖11b)、ASF(圖11c)、AF(圖11d)、FC(圖11e)和常規(guī)化極結(jié)果(圖11i)均得到改善，條帶狀異?；颈幌?，礦區(qū)的正負(fù)磁異常差異較明顯；NT方法(圖11f)的辨識(shí)度增強(qiáng)。整體來(lái)看，正則化改進(jìn)后的七種化極結(jié)果、TT(圖11g)以及HS(圖11h)結(jié)果較為合理，礦區(qū)的化極磁異常特征較明顯，已基本達(dá)到化極目的。

圖10 實(shí)際數(shù)據(jù)九種不同方法的磁異?；瘶O結(jié)果對(duì)比

(a)SF：β=20°，Corr=0.212； (b)DD：d=1，Corr=0.267； (c)ASF：θ′=71°，Corr=0.225； (d)AF：p=1，Corr=0.279； (e)FC：k=2，h=1，Corr=0.283； (f)NT：Q=1，Corr=0.318； (g)TT：M1=5，Corr=0.392； (h)HS：A=2，B=0.1，Corr=0.364； (i)RO：Corr=0.121

圖11 實(shí)際數(shù)據(jù)正則化改進(jìn)后磁異?；瘶O結(jié)果對(duì)比(a)SF； (b)DD； (c)ASF； (d)AF； (e)FC； (f)NT； (g)TT； (h)HS； (i)RO

由于實(shí)際地磁環(huán)境十分復(fù)雜，且更容易受噪聲的影響，因此實(shí)際數(shù)據(jù)處理的難度更大。應(yīng)用改進(jìn)的化極算法能夠在一定程度上降低低緯度磁異常數(shù)據(jù)解釋的難度，同時(shí)，九種不同化極方法能夠在化極結(jié)果上互相驗(yàn)證，這也從另外一個(gè)角度證實(shí)了磁異常化極結(jié)果的可靠性。

6 結(jié)論

(1)基于對(duì)RO的改進(jìn)，局部壓制類(lèi)和全局壓制類(lèi)化極方法能對(duì)不穩(wěn)定的放大區(qū)進(jìn)行高頻壓制，且全局壓制類(lèi)方法的壓制效果要優(yōu)于局部壓制類(lèi)方法。

(2)局部壓制類(lèi)和全局壓制類(lèi)化極方法在實(shí)現(xiàn)對(duì)放大區(qū)高頻壓制的同時(shí)也引入了參數(shù)，并且參數(shù)的取值直接影響化極結(jié)果的精度。通過(guò)計(jì)算化極結(jié)果與NSS的最大相關(guān)系數(shù)可以得到最優(yōu)參數(shù)和最優(yōu)參數(shù)下的化極結(jié)果，提高化極結(jié)果的精度。

(3)參數(shù)最優(yōu)的化極因子在較低緯度仍存在高頻放大效應(yīng)，利用正則化改進(jìn)算法能有效抑制原始信號(hào)中的高頻噪聲并提高化極結(jié)果的信噪比。

本文通過(guò)仿真試驗(yàn)和實(shí)際數(shù)據(jù)實(shí)驗(yàn)對(duì)比分析了九種化極方法的應(yīng)用效果，并給出了相應(yīng)的優(yōu)化方法，最終得到了穩(wěn)定的磁異常化極結(jié)果。在實(shí)際應(yīng)用中，應(yīng)根據(jù)緯度高低、噪聲水平和實(shí)時(shí)化極效果等因素綜合考慮，選取最優(yōu)方法以獲得最優(yōu)化極結(jié)果。

值得注意的是，磁異常數(shù)據(jù)經(jīng)過(guò)化極或化赤處理都能夠降低斜磁化的復(fù)雜性。在低緯度地區(qū)，尤其是磁赤道處，化赤因子是穩(wěn)定的。由于水平磁化的磁異常數(shù)據(jù)解釋難度較大，因此也可以對(duì)低緯度磁異常數(shù)據(jù)進(jìn)行化赤處理，然后倒相180°。

參考文獻(xiàn)

[1] Baranov V.A new method for interpretation of aeromagnetic maps:pseudo-gravimetric anomalies.Geophysics,1957,22(2):359-382.

[2] Silva J B C.Reduction to the pole as an inverse pro-blem and its application to low-latitude anomalies.Geo-physics,1986,51(2):369-382.

[3] Leao J W D,Silva J B C.Discrete linear transformations of potential field data.Geophysics,1989,54(4):497-507.

[4] Bhattacharyya B K.Two-dimensional harmonic analysis as a tool for magnetic interpretation.Geophysics,1965,30(5):829-857.

[5] 柴玉璞.從化極算法誤差方程看各種波數(shù)域低緯度化極方法.石油地球物理勘探,2012,47(3):496-505.

Chai Yupu.Discussion of wavenumber domain RTP methods at low latitudes.OGP,2012,47(3):496-505.

[6] 姚長(zhǎng)利，管志寧，高德章等.低緯度磁異?；瘶O方法——壓制因子法.地球物理學(xué)報(bào),2003,46(5):690-696.

Yao Changli,Guan Zhining,Gao Dezhang et al.Reduction to the pole of magnetic anomalies at low latitude with suppression filter.Chinese Journal of Geophysics,2003,46(5):690-696.

[7] Li X.Magnetic reduction-to-the-pole at low latitudes:Observations and considerations.The Leading Edge,2008,27(8):990-1002.

[8] 姚長(zhǎng)利,黃衛(wèi)寧,張聿文等.直接阻尼法低緯度磁異?；瘶O技術(shù).石油地球物理勘探,2004,39(5):600-606.

Yao Changli,Huang Weining,Zhang Yuwen et al.Reduction to the pole of magnetic anomalies at low latitude with direct damping.OGP,2004,39(5):600-606.

[9] 林曉星,王平.一種改進(jìn)的低緯度磁場(chǎng)化極方法——變頻雙向阻尼因子法.地球物理學(xué)報(bào),2012,55(10):3477-3484.

Lin Xiaoxing,Wang Ping.An improved method for reduction to the pole of magnetic field at low latitude -the method of frequency conversion bidirectional damping factor.Chinese Journal of Geophysics,2012,55(10):3477-3484.

[10] Guo L,Shi L,Meng X.The antisymmetric factor me-thod for magnetic reduction to the pole at low latitudes.Journal of Applied Geophysics,2013,92:103-109.

[11] Zhang H,Marangoni Y R,Hu X et al.NTRTP:A new reduction to the pole method at low latitudes via a nonlinear thresholding.Journal of Applied Geophy-sics,2014,111:220-227.

[12] Mendona C A,Silva J B C.A stable truncated series approximation of the reduction-to-the-pole operator.Geophysics,1993,58(8):1084-1090.

[13] 張培琴,趙群友.低磁緯度區(qū)航磁異常變傾角磁方向轉(zhuǎn)換方法.物探化探計(jì)算技術(shù),1996,18(3):206-214.

Zhang Peiqin,Zhao Qunyou.Methods of the magnetic direction transform of aeromagnetic anomalies with differential inclinations in low magnetic latitudes.Computing Techniques for Geophysical and Geoche-mical Exploration,1996,18(3):206-214.

[14] 張振宇,王剛,胡祥云等.張量CSAMT方法及對(duì)比實(shí)驗(yàn).石油地球物理勘探,2017,52(4):869-874.

Zhang Zhenyu,Wang Gang,Hu Xiangyun et al.Tensor CSAMT technical and experiments.OGP,2017,52(4):869-874.

[15] Wilson H.Analysis of the magnetic gradient tensor.Defence Research Establishment Pacific:Canada Technical Memorandum,1985,85:13-47.

[16] Zeng X,Liu D,Li X et al.An improved regularized downward continuation of potential field data.Journal of Applied Geophysics,2014,106:114-118.

[17] Yin G,Zhang Y,Mi S et al.Calculation of the magne-tic gradient tensor from total magnetic anomaly field based on regularized method in frequency domain.Journal of Applied Geophysics,2016,134:44-54.

[18] 曹靜杰.基于廣義高斯分布和非Lp范數(shù)正則化的地震稀疏盲反褶積.石油地球物理勘探,2016,51(3):428-433.

Cao Jingjie.Seismic sparse blind deconvolution based on generalized Gaussian distribution and non-convex Lpnormal regularization.OGP,2016,51(3):428-433.

[19] 楊國(guó)權(quán),丁鵬程,李振春等.應(yīng)用線性程函方程和整形正則化的三維初至波旅行時(shí)層析.石油地球物理勘探,2017,52(2):264-272.

Yang Guoquan,Ding Pengcheng,Li Zhenchun et al.3D first-arrival traveltime tomography with linearized eikonal and shaping regularization.OGP,2017,52(2):264-272.

[20] Hansen P C,O’Leary D P.The use of the L-curve in the regularization of discrete ill-posed problems.SIAM Journal on Scientific Computing,1993,14(6):1487-1503.