自適應(yīng)正則化濾波化極方法研究

何濤, 王皓, 張義蜜, 王萬銀,2,3*, Colin G. Farquharson,J. Kim Welford

1 長安大學地質(zhì)工程與測繪學院, 西安 710054

2 海洋油氣勘探國家工程研究中心, 北京 100028

3 中國科學院海洋地質(zhì)與環(huán)境重點實驗室, 山東青島 266071

4 Department of Earth Sciences, Memorial University of Newfoundland, St. John′s, AlB 3X5, Canada

0 引言

磁力異常受到斜磁化的影響,對正確認識磁性體造成了困難.所以使用化極處理消除斜磁化的影響,成為利用磁力異常進行解釋的首要環(huán)節(jié).相比于空間域的褶積方法(Baranov,1957),在波數(shù)域(頻率域)進行化極僅是簡單的乘積運算(Bhattacharyya,1965),其表達式簡單且效率更高,但在低緯度地區(qū)會出現(xiàn)不穩(wěn)定性問題.因此眾多學者針對低緯度地區(qū)化極的不穩(wěn)定性問題做了大量的研究,提出了許多方法.在空間域,主要為等效源方法(Silva,1986;郭志宏,1997;駱遙和薛典軍,2009;Li et al.,2014),該方法計算精度高,但是計算效率低.在波數(shù)域(頻率域),主要采用濾波方法(Macleod et al.,1993;姚長利等,2003;Li,2008;石磊等,2012;Guo et al.,2013;Zhang et al.,2014;Stewart,2019;Ribeiro,2020)來提高化極的穩(wěn)定性,其中比較典型的一類就是基于正則化思想的正則化濾波方法,包括維納濾波法(Hansen and Pawlowski,1989;Keating and Zerbo,1996)、雙曲正(余)弦法(張培琴和趙群友,1996)、直接阻尼法(姚長利等,2004;荊磊等,2017)、變頻雙向阻尼法(林曉星和王平,2012)和正則化法(曾小牛等,2016;張琪等,2018).正則化濾波方法是在化極因子中引入正則化濾波函數(shù),來壓制低緯度地區(qū)化極因子的放大作用.該類方法的關(guān)鍵在于正則化濾波函數(shù)的構(gòu)建,前人一般采用正弦(余弦)函數(shù),以保證其連續(xù)性、有限性,但是其中設(shè)計的控制參數(shù)較多,關(guān)鍵參數(shù)選擇標準缺乏理論研究,從而降低了方法的實用性.因此為提高方法的實用性,當前廣泛利用受磁化方向影響較小的磁場轉(zhuǎn)換量與化極結(jié)果的相關(guān)系數(shù)來確定方法的參數(shù)(Dannemiller and Li,2006;Gerovska et al.,2009;Zhang et al.,2014;Rao et al.,2016;Li et al., 2017;Zhang et al.,2018;Hao et al.,2018).常見的受磁化方向影響較小的磁場轉(zhuǎn)換量主要是解析信號振幅(Nabighian,1972)、磁場模量及其相關(guān)轉(zhuǎn)換量(Stavrev and Gerovska,2000;Gerovska and Araúzo-Bravo,2006)以及歸一化磁源強度(Wynn et al.,1975).其中,解析信號振幅較易獲得,應(yīng)用最為廣泛,但是三度體的解析信號振幅仍然會受到磁化方向的強烈影響,此外還受地質(zhì)體的埋深影響(Li,2006;王萬銀,2012).而磁場模量及其相關(guān)轉(zhuǎn)換量和歸一化磁源強度需已知磁場三分量才能有效地降低磁化方向的影響,但是在實際工作中,較多的情況為僅已知磁場總場,在低緯度地區(qū),利用磁場總場求取磁場三分量同樣存在不穩(wěn)定性問題.故通過相關(guān)性確定化極方法的參數(shù)存在局限性,不能從根本上解決參數(shù)選擇的問題.

此外,正則化濾波方法在壓制放大作用的同時,由于其濾波尺度難以把握,所以普遍存在濾波過剩現(xiàn)象,損失了有用信號,降低了該類方法化極結(jié)果的精度.而迭代方法(Strakhov and Devitsyn,1965;徐世浙,2006)是解決位場數(shù)據(jù)處理與轉(zhuǎn)換的不穩(wěn)定性問題的另一重要方法,在低緯度地區(qū)的化極處理中也得到了廣泛應(yīng)用(Li,2008;駱遙和薛典軍,2010;姚長利等,2012;Hao et al.,2018;He et al.,2022),該類方法僅需近似的化極因子,然后利用原始磁力異常進行約束,通過多次計算提高化極結(jié)果的精度.因此將構(gòu)建的正則化化極因子作為一個近似因子,代入迭代過程,可以有效地提高正則化濾波方法的精度.上述所有研究者均聚焦于低緯度地區(qū)具體化極方法的研究,但是關(guān)于其應(yīng)用條件“低緯度地區(qū)”這個概念的界定卻十分模糊,即缺少高、低緯度地區(qū)的分界問題的研究.

綜上,明確低緯度地區(qū)的范圍和研究實用性強、精度高的低緯度地區(qū)化極方法具有重要意義.為此,本文依據(jù)理論研究確定了化極因子振幅引起放大作用的臨界值,進而確定高、低緯度地區(qū)磁化傾角的分界值;其中,在低緯度地區(qū)依據(jù)正則化思想,提出自適應(yīng)正則化濾波化極方法(The adaptive regularization filtering method of reduction-to-the-pole,ARF-RTP),并通過理論研究確定不同磁化傾角條件下該方法的參數(shù)值,提高方法的實用性.同時,結(jié)合迭代法解決正則化濾波作用導(dǎo)致的有用信號損失問題,提高方法的精度,有效地解決低緯度地區(qū)化極的不穩(wěn)定性問題.最后,通過理論模型測試和實際資料處理驗證了該方法的正確性、穩(wěn)定性和實用性.

1 化極因子特征分析

波數(shù)域(頻率域)磁場分量和磁化方向轉(zhuǎn)換因子H(u,v)為(Bhattacharyya,1965)

H(u,v)=

(1)

將(1)式轉(zhuǎn)換至極坐標,得到極坐標下轉(zhuǎn)換因子表現(xiàn)形式H(θ)(Mendon?a and Silva,1993)：

(2)

式中,θ=arctan(v/u).如果僅進行磁化方向轉(zhuǎn)換,式中保持分量方向在轉(zhuǎn)換前后不變即可.同理,磁場分量方向轉(zhuǎn)換也是相同.

當磁場分量方向和磁化方向一致時,(2)式可以簡化為(Mendon?a and Silva,1993)

(3)

磁場分量和磁化方向轉(zhuǎn)換中最常使用的是化極(Reduction-to-the-pole,RTP),即轉(zhuǎn)換到垂直磁化的垂直分量(I2=90°),其轉(zhuǎn)換因子(Mendon?a and Silva,1993)為

(4)

將其分解為振幅ARTP(θ)和相位φRTP(θ)形式(荊磊等,2017)：

(5)

在波數(shù)域(頻率域)轉(zhuǎn)換時,振幅相乘,而相位則是執(zhí)行相加運算.(5)式中相位滿足-180°≤φRTP(θ)≤180°,在轉(zhuǎn)換中僅起到改變原始波譜相位的作用,不會引起不穩(wěn)定性問題,而振幅ARTP(θ)值域與磁化傾角I1相關(guān),其參與的相乘運算會導(dǎo)致不穩(wěn)定性問題.為了進一步說明化極因子振幅和相位特征,做出其相應(yīng)的等值線圖(圖1).

圖1 全傾角化極因子的振幅與相位等值線圖(a) D=0°振幅; (b) D=0°相位; (c) D=45°振幅; (d) D=45°相位. 振幅圖中等值線間隔按照對數(shù)函數(shù)(ln函數(shù))值分布規(guī)律進行設(shè)置;振幅圖中最大值為無窮大,為了便于成圖,設(shè)置色標最大值為1000;相位圖中等值線間隔為30.

對比圖1a和圖1c可以看出化極因子振幅在南、北半球變化特征相同,呈現(xiàn)明顯的對稱性,最小值為1;各磁化傾角的振幅均在θ-D=±90°處達到最大值,隨著磁化傾角的減小,最大值逐漸增大,而磁化偏角決定了振幅放大的位置,但它本身不會影響振幅的放大效果.對比圖1b和圖1d可以看出相位的范圍始終在-180°～180°之間,在θ-D=±90°時相位為0,在θ-D=±90°兩側(cè)由正(負)相位變?yōu)樨?正)相位,隨著磁化傾角的減小,相位曲線變得陡峭,而磁化偏角決定了相位反相的位置.

圖1中,各振幅均以2為界(圖中白色實線)呈現(xiàn)不同的特征.當振幅小于2時,等值線呈現(xiàn)非圓形特征,變化較緩慢,隨著幅值逐漸增大,形狀收斂于以放大中心(θ-D=±90°,0°)為圓心的圓形;當振幅等于2時,等值線形成以(θ-D=±90°,0°)為圓心,半徑為45°的圓形;當振幅大于2時,等值線依舊保持以(θ-D=±90°,0°)為圓心的圓形特征,但是其等值線代表的值急劇增大,至圓心趨于無窮(各振幅圖中最大值為無窮大,但網(wǎng)格化后顯示為特征值,所以振幅圖的色標顯示有界,在這里為了便于成圖,最大值設(shè)置為1000),說明在幅值2形成的圓內(nèi),無論是縱軸I還是橫軸θ,其趨近于圓心(θ-D=±90°,0°)方向的相同變化量所引起振幅的變化量也是相同的,這種變化接近圓心時趨于無窮.

為了更清晰的研究化極因子振幅變化的特征,在這里使用總水平導(dǎo)數(shù)(Cordell,1979)研究其變化率,化極因子振幅總水平導(dǎo)數(shù)如圖2所示.

圖2顯示振幅總水平導(dǎo)數(shù)以振幅2為界(圖中白色線條).當振幅小于2時,振幅無論是在橫軸I方向還是縱軸θ方向變化率均較小,其總水平導(dǎo)數(shù)值小于0.1;而當振幅大于2時,振幅變化劇烈,總水平導(dǎo)數(shù)值大于0.1且急劇增大,反映了振幅正向增大的變化趨勢,在(θ-D=±90°,0°)處振幅變化率同振幅一樣趨于無窮.

化極因子振幅無界,而最小值僅為1,幅值的量級變化非常大,討論中受大值的影響,弱小信號會被湮沒.而且在總水平導(dǎo)數(shù)計算過程中,其計算結(jié)果受網(wǎng)格間距影響較大,為使結(jié)果更具普適性,我們希望得到總水平導(dǎo)數(shù)特征的分界性而非幅值的分界,所以結(jié)合化極的反變換過程進行討論,即由化極場轉(zhuǎn)換到任意低緯度磁場的過程,該轉(zhuǎn)換因子的振幅ARTP→Lowfield(θ)為

ARTP→Lowfield(θ)=cos2(I1)·cos2(θ-D1)+sin2(I1),

(6)

該因子與ARTP(θ)互為倒數(shù),幅值在0～1之間,同樣對其進行求導(dǎo)分析.

圖2 全傾角化極因子振幅的總水平導(dǎo)數(shù)等值線圖(a) D=0°振幅總水平導(dǎo)數(shù); (b) D=45°振幅總水平導(dǎo)數(shù). 等值線間隔按照對數(shù)函數(shù)(ln函數(shù))值分布規(guī)律進行設(shè)置;圖中白色線條表示振幅值為2的等值線.

圖3表示ARTP→Lowfield(θ)的振幅以及其總水平導(dǎo)數(shù)圖.該振幅(圖3a,圖3c)和化極因子振幅特征相同,在南、北半球呈現(xiàn)明顯的對稱性,最小值在(θ-D=±90°,0°)處,幅值為0.分別求其總水平導(dǎo)數(shù)得到圖3b,圖3d,其最大值反映了ARTP→Lowfield(θ)最大變化率,恰在ARTP→Lowfield(θ)=0.5(圖中白色實線)處,當ARTP→Lowfield(θ)>0.5時,隨著ARTP→Lowfield(θ)值減小,總水平導(dǎo)數(shù)值增大,說明其變化率增大,對應(yīng)的ARTP(θ)<2;當ARTP→Lowfield(θ)=0.5時,其對應(yīng)的ARTP(θ)=2;當ARTP→Lowfield(θ)<0.5時,隨著ARTP→Lowfield(θ)值減小,總水平導(dǎo)數(shù)值也減小,說明其變化率減小,其對應(yīng)的ARTP(θ)>2.所以0.5是ARTP→Lowfield(θ)幅值變化特征的分界,與之對應(yīng)的2是ARTP(θ)幅值變化特征的分界.

結(jié)合ARTP(θ)的特征、ARTP(θ)的變化率特征、ARTP→Lowfield(θ)的特征以及ARTP→Lowfield(θ)的變化率特征,可以看出當ARTP(θ)≤2時,雖然其值逐漸變大但是其變化率較小,振幅最大值與最小值的偏差在0～1之間,保持穩(wěn)定狀態(tài);當ARTP(θ)>2時,ARTP(θ)其值依舊保持變大趨勢,但是其變化率較大,振幅在θ-D=±90°左右兩側(cè)對稱的扇形范圍內(nèi)開始大于2,而且隨著|I|的減小,這個扇形范圍增大,扇形范圍內(nèi)振幅的值也增大,振幅最大值與最小值的偏差大于1,并趨于無窮,為不穩(wěn)定狀態(tài).綜上,幅值2為化極因子振幅穩(wěn)定與非穩(wěn)定狀態(tài)的分界,其對應(yīng)的磁化傾角|I|=45°,而|I|<45°的范圍則為化極過程中存在不穩(wěn)定問題的低緯度地區(qū).

圖3 全傾角化極場轉(zhuǎn)換低緯度地區(qū)磁場轉(zhuǎn)換因子的振幅與總水平導(dǎo)數(shù)等值線圖(a) D=0°振幅; (b) D=0°振幅總水平導(dǎo)數(shù); (c) D=45°振幅; (d) D=45°振幅總水平導(dǎo)數(shù). 總水平導(dǎo)數(shù)圖中白色線條表示振幅值為0.5的等值線;振幅圖中等值線間隔為0.1;水平導(dǎo)數(shù)圖中等值線間隔為0.001.

2 自適應(yīng)正則化濾波化極方法

位場數(shù)據(jù)的處理和轉(zhuǎn)換通常存在這樣一類問題,其轉(zhuǎn)換因子是如下的分數(shù)形式：

(7)

(8)

其中,對正則化濾波函數(shù)R(θ)的基本要求是該函數(shù)的引入僅是提高轉(zhuǎn)換因子的穩(wěn)定性,但是不能改變轉(zhuǎn)換因子的基本特征.那么簡單、高效的正則化濾波函數(shù)的求取是提高正則化方法實用性的關(guān)鍵.

2.1 自適應(yīng)正則化濾波方法

化極因子的振幅形式類同于公式(7),其放大作用是由振幅的分母趨于零所引起的,因此可以利用正則化思想,通過在振幅的分母上引入正則化濾波函數(shù),去壓制放大作用.

在這里我們利用化極因子振幅構(gòu)建正則化濾波函數(shù)：

(9)

式中,α為正則化因子,是一個非負實數(shù),用來調(diào)節(jié)正則化濾波函數(shù)作用的大小.將(9)式引入到化極因子振幅的分母并化簡,得到正則化的化極因子振幅：

Regulation_ARTP(θ)=

(10)

此時構(gòu)建的正則化的化極因子振幅的分子與分母變化規(guī)律保持一致,可以實現(xiàn)自適應(yīng)抑制化極因子振幅的放大作用,而且無須設(shè)置多余參數(shù),在使用中僅需確定正則化因子α的值即可.

前文分析了化極因子的特征,其放大作用是由于振幅大于2所引起的,而振幅在1～2之間是穩(wěn)定的,所以為了盡可能地保留化極因子本身的有效信號,同時提高穩(wěn)定性,我們采用區(qū)域正則化,此時構(gòu)建的正則化振幅AARF-RTP(θ)為如下形式：

(11)

以上就是自適應(yīng)正則化濾波化極方法的基本原理,其化極因子為如下形式：

ARF_HRTP(θ)=AARF-RTP(θ)·eiφRTP(θ),

(12)

其中,AARF-RTP(θ)表達見式(11),φRTP(θ)表達見式(5).因此,利用自適應(yīng)正則化濾波化極方法計算化極場的具體公式如下：

RTP(θ)=[AARF-RTP(θ)·AΔT(θ)]·ei[φRTP(θ)+φΔT(θ)],

(13)

其中,RTP(θ)表示化極場波譜;AΔT(θ)為原始磁場振幅,與上文構(gòu)建的自適應(yīng)正則化振幅相乘;φΔT(θ)為原始磁場相位,與化極因子本身的相位相加.

2.2 正則化因子α的確定

本文通過理論研究確定正則化因子α的值.當α取不同值,可以得到不同θ處的AARF-RTP(θ)值,選擇其最大值構(gòu)建如下函數(shù)：

f(α)=max[AARF-RTP(θ)].

(14)

前文分析可知,當振幅大于2時,存在放大作用,因此以2作為約束,使f(α)函數(shù)滿足下列條件：

f(α)≤2.

(15)

滿足上述條件的α存在無窮多個,但是α越大,濾波作用越強,可能會出現(xiàn)過濾波問題,因此選擇該條件下的最小α值,既保證了穩(wěn)定轉(zhuǎn)換,同時減小過濾波的影響.根據(jù)上述分析計算得到不同磁化傾角對應(yīng)的最佳的α值(見表1).

表 1 各磁化傾角對應(yīng)的α值表Table 1 α corresponding to each magnetization inclination

表1給出了各整數(shù)磁化傾角對應(yīng)的α值,當磁化傾角為小數(shù)時,選擇其相鄰的兩個整數(shù)磁化傾角對應(yīng)的α值中的最大值即可,以保證完全壓制.

將各傾角確定的α值代入式(11),得到全傾角自適應(yīng)正則化濾波化極因子的振幅特征圖(圖4,可以看出該正則化振幅保持了其原來振幅對稱性和連續(xù)型的特征,而且當|I|<45°時,正則化振幅最大值為2(圖4a),在θ-D=±90°處達到極小值,向兩側(cè)對稱的增大,直至達到最大值2之后又開始減小(圖4b),但是部分傾角的最小值小于1,說明該正則化振幅可實現(xiàn)穩(wěn)定轉(zhuǎn)換,但是會引起信號損失.

圖4 全傾角自適應(yīng)正則化濾波化極因子振幅圖(a) 振幅(D=0°); (b) 不同磁化傾角振幅對比(D=0°). 等值線間隔為0.1.

2.3 迭代法提高化極精度

濾波類化極方法的化極結(jié)果適用于定性解釋,不利于定量解釋(柴玉璞和萬海珍,2020),其主要原因是在提高化極穩(wěn)定性的同時,導(dǎo)致化極損失效應(yīng)增強,進而影響化極結(jié)果的精度.而本次研究的自適應(yīng)正則化濾波化極方法也屬于濾波類方法,雖然其盡可能地保留了化極因子本身的有效信號,但依舊會產(chǎn)生損失效應(yīng).所以結(jié)合迭代法,以原始磁力異常作為約束,通過多次計算削弱損失效應(yīng),提高化極結(jié)果精度.

位場中數(shù)據(jù)轉(zhuǎn)換的迭代法收斂通式為(姚長利等,2012)

|1-φ(θ)κ-1(θ)|<1,

(16)

其中,φ(θ)為近似轉(zhuǎn)換因子,在這里為式(11)表示的正則化振幅因子;κ(θ)為真實轉(zhuǎn)換因子,即式(5)中的真實振幅因子,將其代入(16)式得到：

(17)

式中0

具體迭代流程(圖5)如下：

圖5 迭代法化極流程圖

(1)根據(jù)式(13),進行自適應(yīng)正則化濾波化極,得到初始化極磁力異常(RTP(1)).

(2)利用化極磁力異常(RTP(1))計算磁力異常(ΔT(1)).

(3)利用原始磁力異常ΔT與計算磁力異常ΔT(1)求取殘差,即δΔT(1)=ΔT-ΔT(1).

(4)判斷殘差δΔT(1)是否滿足精度控制條件：δΔT(1)<ε(ε為迭代控制誤差).滿足精度條件,輸出化極場;不滿足精度條件,返回到步驟(1)對殘差δΔT(1)進行正則化濾波化極處理,然后利用殘差化極場δRTP(1)對初始化極磁力異常RTP(1)進行校正,得到RTP(2)=RTP(1)+δRTP(1).

依次重復(fù)以上步驟直到殘差δΔT(k)滿足控制要求輸出化極異常RTP(k)即可.

化極轉(zhuǎn)換是磁化方向和磁場分量同時進行轉(zhuǎn)換的問題,因此當我們僅執(zhí)行分量轉(zhuǎn)換或者磁化方向單一轉(zhuǎn)換時,均可按照上述的原理進行處理.

3 模型試算

本研究分別設(shè)計單一模型和復(fù)雜模型模擬實際地質(zhì)體測試方法的正確性和穩(wěn)定性.采用目前較為廣泛應(yīng)用的商業(yè)軟件Geosoft的Oasis montaj軟件(8.4版本)和GeoProbe軟件(4.0版本)與本文ARF-RTP方法進行化極效果對比.Oasis montaj軟件中化極方法為偽傾角法;GeoProbe軟件中化極方法選擇為雙曲余弦濾波法,上述兩個軟件所使用的相關(guān)參數(shù)為多次測試后化極效果較為理想的參數(shù).通過求取理論化極磁力異常與計算得到的化極磁力異常之間的均方根誤差：

(18)

相對均方根誤差：

(19)

說明各個方法的精度.其中,M,N為點數(shù)和線數(shù),RTP(i,j)表示第(i,j)點上正演的理論化極磁力異常值,RTPcal(i,j)表示第(i,j)點上計算得到的化極磁力異常值,RTPmax表示理論化極磁力異常的最大值,RTPmin表示理論化極磁力異常的最小值.

下文各圖中白色線框表示形體的平面位置,文中計算得到的化極磁力異常與其對應(yīng)的理論化極磁力異常的色標以及等值線間隔均保持一致.

3.1 單體模型測試

單體模型選擇由Hansen和Pawlowski(1989)設(shè)計的模型,該模型被眾多學者用來測試各方法的化極效果,模型具體參數(shù)見表2(x1、x2,y1、y2,z1、z2分別為各直立六面體沿x,y,z方向的范圍),設(shè)置觀測界面為z=0.0 m(z坐標方向鉛垂向下為正).

圖6a為單體模型正演的I=0°磁化磁力異常;圖6b為正演的I=90°磁化磁力異常,異常值在-22.43～62.40 nT之間,異常幅值為84.83 nT.為了更加直觀的對比各個方法的化極效果,設(shè)置了通過主異常體的A-A′剖面(x=31.5 m)與B-B′剖面(y=31.5 m)(圖6b中紅色線條)對比化極結(jié)果.

表2 單體模型參數(shù)及測網(wǎng)參數(shù)表Table 2 Parameters and measurement network of the simple model

圖6 單體模型正演磁力異常圖(a) I=0°磁化磁力異常; (b) I=90°磁化磁力異常. 圖中等值線間隔均為5 nT.

圖7 單體模型磁力異?；瘶O結(jié)果圖(a) Oasis montaj軟件化極結(jié)果; (b) GeoProbe軟件化極結(jié)果; (c) ARF-RTP法化極結(jié)果; (d) A-A′剖面化極結(jié)果對比圖; (e) B-B′剖面化極結(jié)果對比圖. 圖中白色線框表示形體平面位置,等值線間隔均為5 nT.

圖7為圖6a表示的磁力異常的化極結(jié)果.Oasis montaj軟件中Amplitude correction inclination選擇為20°,GeoProbe軟件中濾波參數(shù)AS選擇為12,BS選擇為0.0001,本文ARF-RTP法中α值選擇0.242.

Oasis montaj軟件(圖7a)和GeoProbe軟件(圖7b)的化極結(jié)果在異常體正上方與理論化極異常(圖6b)特征相差較大,受沿磁化偏角方向的拉伸現(xiàn)象的影響,并未完全形成與理論化極異常一致的圓弧形圈閉異常特征.本文ARF-RTP法計算的化極結(jié)果(圖7c)明顯減弱了沿磁化偏角方向的拉伸現(xiàn)象的影響,基本恢復(fù)了異常體正上方的異常特征.通過A-A′剖面(圖7d)和B-B′剖面(圖7e)整體可以看出,本文ARF-RTP法化極結(jié)果更接近理論化極異常,其次為GeoProbe軟件,Oasis montaj軟件效果最差.Oasis montaj軟件、GeoProbe軟件以及ARF-RTP法化極結(jié)果的均方誤差和相對均方誤差統(tǒng)計見表3,ARF-RTP法精度明顯優(yōu)于上述兩個軟件.細節(jié)方面,在A-A′剖面(圖7d)方向ARF-RTP法、GeoProbe軟件、Oasis montaj軟件化極結(jié)果形態(tài)與理論異常形態(tài)基本一致,但是Oasis montaj軟件化極結(jié)果與理論化極磁力異常存在背景差,而ARF-RTP法與GeoProbe軟件化極結(jié)果沒有表現(xiàn)出明顯的背景差;在B-B′剖面(圖7e)方向ARF-RTP法化極結(jié)果與理論異常形態(tài)基本一致,GeoProbe軟件、Oasis montaj軟件化極結(jié)果較理論異常光滑,表明GeoProbe軟件和Oasis montaj軟件方法在壓制放大作用的同時,對垂直于磁化偏角方向的有效信息進行了較強的濾波.

表3 單體模型磁力異?；瘶O結(jié)果誤差統(tǒng)計表Table 3 Error statistics of RTP results of magnetic anomaly of the simple model

3.2 復(fù)雜模型測試

復(fù)雜模型選擇由He等(2022)使用的模型,該模型利用五個大小、埋深各不相同的直立六面體構(gòu)成.模型具體參數(shù)見表4,各參數(shù)代表意義與上述單體模型相同,設(shè)置觀測界面為z=0.0 km.

表4 復(fù)雜模型參數(shù)及測網(wǎng)參數(shù)表Table 4 Parameters and measurement network of the complex model

圖8 復(fù)雜模型正演磁力異常圖(a) I=0°磁化磁力異常; (b) I=10°磁化磁力異常; (c) I=30°磁化磁力異常; (d) I=10°磁化含噪聲磁力異常; (e) I=90°磁化磁力異常. 圖中白色線框表示形體平面位置,等值線間隔為50 nT.

圖9 復(fù)雜模型磁力異?；瘶O結(jié)果(a) I=0°磁力異常Oasis montaj軟件化極結(jié)果; (b) I=0°磁力異常GeoProbe軟件化極結(jié)果; (c) I=0°磁力異常ARF-RTP法化極結(jié)果; (d) I=10°磁力異常Oasis montaj軟件化極結(jié)果; (e) I=10°磁力異常GeoProbe軟件化極結(jié)果; (f) I=10°磁力異常ARF-RTP法化極結(jié)果; (g) I=30°磁力異常Oasis montaj軟件化極結(jié)果; (h) I=30°磁力異常GeoProbe軟件化極結(jié)果; (i) I=30°磁力異常ARF-RTP法化極結(jié)果. 圖中白色線框表示形體平面位置,等值線間隔均為50 nT.

圖10 復(fù)雜模型磁力異常不同方法化極結(jié)果對比(a) A-A′剖面I=0°磁力異?；瘶O結(jié)果對比圖; (b) B-B′剖面I=0°磁力異?；瘶O結(jié)果對比圖; (c) A-A′剖面I=10°磁力異常化極結(jié)果對比圖; (d) B-B′剖面I=10°磁力異?；瘶O結(jié)果對比圖; (e) A-A′剖面I=30°磁力異常化極結(jié)果對比圖; (f) B-B′剖面I=30°磁力異?；瘶O結(jié)果對比圖.

圖11 復(fù)雜模型含噪聲磁力異?；瘶O結(jié)果圖(a) Oasis montaj軟件化極結(jié)果; (b) GeoProbe軟件化極結(jié)果; (c) ARF-RTP法化極結(jié)果; (d) A-A′剖面化極結(jié)果對比圖; (e) B-B′剖面化極結(jié)果對比圖. 文中所加噪聲選擇均值為0 nT,標準偏差為3 nT的高斯噪聲,圖中白色線框表示形體平面位置,等值線間隔均為50 nT.

圖8(a,b,c)分別表示復(fù)雜模型正演的I=0°,I=10°,I=30°磁化磁力異常圖;圖8d表示對圖8b表示的磁力異常加均值為0 nT,標準偏差為3 nT高斯噪聲的含噪聲磁力異常;圖8e為正演的I=90°磁化磁力異常,異常值在-144.11-1092.29 nT之間,異常幅值為1236.40 nT.選擇A-A′剖面(x=7 km)與B-B′剖面(y=12 km)(圖8e中紅色線條)進行化極結(jié)果剖面對比.

圖9為圖8a、圖8b、圖8c表示的磁力異常的化極結(jié)果.其中,在I=0°時,Oasis montaj軟件中Amplitude correction inclination選擇為10°,GeoProbe軟件中濾波參數(shù)AS選擇為6,BS選擇為0.001,本文ARF-RTP方法中α值選擇0.242;在I=10°時,Oasis montaj軟件中Amplitude correction inclination選擇為15°,GeoProbe軟件中濾波參數(shù)AS選擇為12,BS選擇為0.001,本文ARF-RTP方法中α值選擇0.250;在I=30°時,Oasis montaj軟件中Amplitude correction inclination選擇為30°,GeoProbe軟件中濾波參數(shù)AS選擇為12,BS選擇為0.0001,本文ARF-RTP方法中α值選擇0.249.

當I=0°與I=10°時,Oasis montaj軟件(圖9a,圖9d)和GeoProbe軟件(圖9b,圖9e)的化極結(jié)果中主體異常特征基本恢復(fù),但是在形體的南北兩側(cè)均存在沿磁化偏角方向的拉伸現(xiàn)象,形成條帶狀虛假異常.而ARF-RTP法計算的化極結(jié)果(圖9c)中沿磁化偏角方向的拉伸現(xiàn)象明顯減弱,在I=10°時,其化極結(jié)果(圖9f)已不存在拉伸現(xiàn)象,異常形態(tài)更接近理論異常形態(tài).當I=30°時,Oasis montaj軟件(圖9g)與GeoProbe軟件(圖9h)化極結(jié)果中已無明顯的沿磁偏角方向的拉伸現(xiàn)象,與理論異常形態(tài)較為接近,ARF-RTP法化極結(jié)果(圖9i)與理論異常在形態(tài)和幅值上基本保持一致.不同磁化傾角條件下,Oasis montaj軟件、GeoProbe軟件以及ARF-RTP法的均方誤差和相對均方誤差統(tǒng)計見表5,ARF-RTP法精度明顯優(yōu)于上述兩個軟件,而且磁化傾角越小,這種優(yōu)勢越明顯.

細節(jié)方面,與單體模型相同,在A-A′剖面(圖10a,圖10c)方向ARF-RTP法、GeoProbe軟件、Oasis montaj軟件化極結(jié)果形態(tài)與理論異常形態(tài)基本一致,但是在異常體上方(10 km至22 km之間)Oasis montaj軟件與GeoProbe軟件化極結(jié)果的幅值整體偏小,似乎存在背景差的影響,而ARF-RTP法化極結(jié)果與理論異常曲線基本重合;在B-B′剖面(圖10b,圖10d)方向ARF-RTP法化極結(jié)果與理論異常形態(tài)基本一致,GeoProbe軟件、Oasis montaj軟件化極結(jié)果較理論異常光滑,尤其在異常局部凸起和凹陷處(3～10 km之間)這種光滑現(xiàn)象尤為明顯,反映了GeoProbe軟件和Oasis montaj軟件方法在壓制垂直于磁化偏角方向的放大作用的同時,會對有效信息進行濾波.當I=30°時,在A-A′剖面(圖10e)與B-B′剖面(圖10f)處,ARF-RTP法、GeoProbe軟件、Oasis montaj軟件化極結(jié)果曲線與理論異常曲線基本重合.

圖11為圖8d所示的含噪聲磁力異常的化極結(jié)果.Oasis montaj軟件中Amplitude correction inclination選擇為15°,GeoProbe軟件中濾波參數(shù)AS選擇為12,BS選擇為0.001,本文ARF-RTP方法中α值選擇0.250.添加噪聲后,Oasis montaj軟件(圖11a)、GeoProbe軟件(11b)和ARF-RTP法(圖11c)計算的化極結(jié)果的異常形態(tài)與不含噪聲磁力異常化極結(jié)果基本是一致的,僅是受噪聲影響異常曲線存在鋸齒狀波動.剖面圖11d和圖11e顯示存在噪聲時,ARF-RTP法的化極結(jié)果仍然更接近理論化極場,而且異常曲線波動較小.各方法均方誤差與相對均方誤差統(tǒng)計見表5,當磁力異常含噪聲時,ARF-RTP法對噪聲可以進行較好的壓制,化極結(jié)果的精度依舊明顯優(yōu)于Oasis montaj軟件和GeoProbe軟件.

表5 復(fù)雜模型磁力異?；瘶O結(jié)果誤差統(tǒng)計表Table 5 Error statistics of RTP results of magnetic anomaly of the complex model

4 實際資料處理

為了測試方法在不同低緯度地區(qū)的實際資料處理中的實用性,選擇位于磁赤道附近的Magur Islands研究區(qū)ΔT磁力異常(磁化傾角為2.78°,磁化偏角為3.66°)和緯度相對較高的惠北凹陷研究區(qū)ΔT磁力異常(磁化傾角29.0°,磁化偏角為-1.43°)分別進行化極處理.

Magur Islands研究區(qū)位于卡羅琳板塊與太平洋板塊交界處,圖12為其ΔT磁力異常圖,該數(shù)據(jù)來源于EMAG2(V2)地磁網(wǎng)格(Earth Magnetic Anomaly Grid,2-arc-minute resolution).Magur Islands是卡羅琳板塊熱點的產(chǎn)物,屬于火山巖島嶼,具有高磁特性(Wu et al., 2016),但是ΔT磁力異常呈現(xiàn)明顯的低磁異常特征,而高磁異常主要分布在火山南北兩側(cè),這是典型的磁赤道附近磁性體磁力異常的表現(xiàn),這種異常特征與已知地質(zhì)認識不符,需對其進行化極處理.

圖12 Magur Islands研究區(qū)ΔT磁力異常圖圖中白色線條表示Magur Islands平面位置,等值線間隔為20 nT.

對圖12所表示ΔT磁力異常進行化極處理得到化極磁力異常(圖13).其中,Oasis montaj軟件中Amplitude correction Inclination選擇為60°,GeoProbe軟件中濾波參數(shù)AS選擇為2,BS選擇為0.1,本文ARF-RTP方法中α值選擇0.242.Oasis montaj軟件(圖13a)、GeoProbe軟件(圖13b)和ARF-RTP法(圖13c)三類方法化極結(jié)果的異常特征大致相似,在Magur Islands處均對應(yīng)高磁異常,整體呈現(xiàn)NE向高低值相間分布的條帶特征,是太平洋洋殼磁條帶的典型表現(xiàn),與已知地質(zhì)認識相符合.為了更加清楚的表現(xiàn)其化極結(jié)果是否穩(wěn)定,即是否存在近NS向虛假異常,對上述化極結(jié)果進行分離,得到分離后的剩余化極磁力異常.Oasis montaj軟件(圖13d)與GeoProbe軟件(圖13e)的剩余化極磁力異常在研究區(qū)西部以及北部存在明顯的近NS向的條帶狀異常,這種近NS向異常并非構(gòu)造特征的表現(xiàn),其主要是由低緯度地區(qū)化極的不穩(wěn)定所引起的,由于近NS向的條帶狀虛假異常的影響,使分離的剩余化極場表現(xiàn)為NNE向與NNW向,而非研究區(qū)地質(zhì)構(gòu)造的NE與NW向.而ARF-RTP法的剩余化極磁力異常(圖13f)并未出現(xiàn)近NS向的條帶狀異常,走向為NE向與NW向,與已知地質(zhì)構(gòu)造認識相符.剩余化極磁力異常的特征說明Oasis montaj軟件與GeoProbe軟件穩(wěn)定性較差,而本文ARF-RTP方法穩(wěn)定性更高.

惠北凹陷研究區(qū)位于我國著名的含油氣盆地珠江口盆地內(nèi)部,其ΔT磁力異常如圖14所示(圖中白色實線表示研究區(qū)內(nèi)構(gòu)造單元界線),該數(shù)據(jù)來源于南海北部海域的航測ΔT磁力異常平面圖.

對圖14所表示ΔT磁力異常進行化極處理得到化極磁力異常(圖15).其中,Oasis montaj軟件中Amplitude correction Inclination選擇為35°,GeoProbe軟件中濾波參數(shù)AS選擇為12,BS選擇為0.005,本文ARF-RTP方法中α值選擇0.25.化極之后,異常正值向北偏移,GeoProbe軟件(圖15b)和ARF-RTP法(圖15c)化極結(jié)果的異常特征大致相似,但是Oasis montaj軟件化極結(jié)果(圖15a)與二者相差較大,在惠州凹陷南部仍然是低磁異常,而GeoProbe軟件和ARF-RTP方法在該處已過渡至高磁異常,可能是Oasis montaj軟件中在壓制放大作用的過程中使用了大于實際磁化傾角的Amplitude correction Inclination導(dǎo)致異常正值向北偏移不足.同樣為了檢驗其化極結(jié)果是否穩(wěn)定,對化極磁力異常分別進行分離,得到分離后的剩余化極磁力異常.Oasis montaj軟件(圖15d)、GeoProbe軟件(圖15e)與ARF-RTP法(圖15f)化極結(jié)果分離得到的剩余化極磁力異常特征相近,均未表現(xiàn)出近NS向的條帶狀特征,說明三類方法的化極結(jié)果均是穩(wěn)定的.

5 結(jié)論及建議

本文通過理論研究得到,在化極處理過程中化極因子振幅穩(wěn)定與否的臨界值為2,與之對應(yīng)的磁化傾角為±45°;當振幅大于2時,磁化傾角在-45°～45°之間,該范圍為低緯度區(qū)域,化極處理存在不穩(wěn)定性問題.為實現(xiàn)低緯度地區(qū)穩(wěn)定地自適應(yīng)化極處理,本文依據(jù)正則化思想,利用化極因子振幅構(gòu)建正則化濾波函數(shù),進而提出了低緯度地區(qū)自適應(yīng)正則化濾波化極方法(ARF-RTP),在該方法中以臨界值2為約束,經(jīng)過理論推導(dǎo)得到各磁化傾角對應(yīng)的正則化因子的值,提高了方法的實用性;同時結(jié)合迭代通過多次計算提高了方法的精度,有效地解決了低緯度地區(qū)化極的不穩(wěn)定性問題.通過理論模型試驗和實際資料處理驗證了本文提出的低緯度自適應(yīng)正則化濾波化極方法的正確性、穩(wěn)定性和實用性,而且其處理效果明顯優(yōu)于目前主流重、磁處理軟件(GeoProbe和Oasis montaj).

本次研究是以固定磁化方向條件為前提,后續(xù)還需研究該方法在變磁化方向條件下的應(yīng)用.

圖13 Magur Islands研究區(qū)化極磁力異常圖(a) Oasis montaj軟件化極結(jié)果; (b) GeoProbe軟件化極結(jié)結(jié)果; (c) ARF-RTP方法化極結(jié)果; (d) 圖13a的剩余化極磁力異常; (e) 圖13b的剩余化極磁力異常; (f) 圖13c的剩余化極磁力異常. 圖中白色線條表示Magur Islands平面位置,化極磁力異常等值線間隔均為30 nT,剩余化極磁力異常等值線間隔為5 nT.

圖14 惠州凹陷研究區(qū)ΔT磁力異常圖圖中白色線條表示構(gòu)造單元界線,等值線間隔為10 nT.

致謝感謝論文評審專家提出的修改意見,也感謝本文編輯對論文的加工和修改.