區(qū)域航磁異常匯編數(shù)據(jù)的長(zhǎng)波成分提取、評(píng)價(jià)與替換技術(shù)

張攀， 杜勁松,2,3*， 王震， 楊明豫， 陳超,2

1 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢 430074 2 地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074 3 地質(zhì)過程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室， 武漢 430074

0 引言

地磁場(chǎng)是地球固有的天然物理場(chǎng)之一，其在固體地球科學(xué)、海洋科學(xué)、環(huán)境科學(xué)、資源勘探、工程勘察以及軍事與國(guó)防等諸多領(lǐng)域發(fā)揮著重要作用(管志寧, 2005; Hinze et al., 2013; 楊華和梁月明, 2013; 熊盛青等, 2014, 2016; Li and Wang, 2018).在眾多的地磁場(chǎng)測(cè)量方式中，航空磁測(cè)相比地面磁測(cè)具有較高的測(cè)量效率，且不受水域、森林、沼澤、沙漠和高山等特殊地貌和復(fù)雜地形的限制，同時(shí)由于飛行是在距地表一定的高度(如50 m～5 km)進(jìn)行的，從而減弱了地表磁性不均勻體與環(huán)境干擾等影響(熊盛青, 2020)；此外，相對(duì)衛(wèi)星磁測(cè)(Zhou et al., 2018)而言，航空磁測(cè)由于距離巖石圈較近，對(duì)于巖石圈磁場(chǎng)的中、短波長(zhǎng)成分具有較高的敏感度與測(cè)量精度.因此，各個(gè)國(guó)家與地區(qū)的區(qū)域巖石圈磁場(chǎng)調(diào)查一般均采用航空測(cè)量方式(Ravat et al., 2009; 熊盛青等, 2013; Xiong et al., 2016; Golynsky et al., 2018; Djomani et al., 2019).

但是，各個(gè)國(guó)家和地區(qū)編制的航磁異常圖及其數(shù)據(jù)集一般均是由多個(gè)數(shù)據(jù)子集匯編而成，這些子數(shù)據(jù)集之間的測(cè)量時(shí)間不同、儀器與平臺(tái)各異，具有不同的飛行高度、磁場(chǎng)測(cè)量分辨率和精度以及測(cè)點(diǎn)定位精度，而且不同時(shí)期采用的數(shù)據(jù)校正方式以及相關(guān)模型(例如正常場(chǎng)校正所采用的主磁場(chǎng)模型的分辨率和精度各異)和參數(shù)也可能不相同.即使上述問題均不存在，由于地球主磁場(chǎng)隨時(shí)間變化，一方面巖石圈感應(yīng)磁化強(qiáng)度及感應(yīng)磁場(chǎng)也將發(fā)生變化(Jackson, 2007; Thébault et al., 2009; Sebera et al., 2019)，另一方面即使巖石圈磁性結(jié)構(gòu)不發(fā)生任何變化，總磁場(chǎng)強(qiáng)度異常即ΔT磁力異常也會(huì)由于主磁場(chǎng)變化而發(fā)生改變(袁曉雨等, 2015; 孫石達(dá)等, 2020).因此，區(qū)域航磁異常匯編數(shù)據(jù)往往存在低頻或長(zhǎng)波成分不可靠性(傅敖云, 1984; Ravat et al., 2003).這雖然對(duì)于巖石圈淺部的、局部的磁性體與地質(zhì)構(gòu)造等研究的影響較小，但是對(duì)于理解深部的、區(qū)域性巖石圈磁性構(gòu)造格架及其大地構(gòu)造與動(dòng)力學(xué)意義具有較大影響.

相比之下，衛(wèi)星磁測(cè)可以在較短的時(shí)間內(nèi)比較均勻地覆蓋全球，對(duì)地磁場(chǎng)的測(cè)量精度與分辨率也在逐步提升，隨著數(shù)據(jù)處理與建模技術(shù)的快速發(fā)展，目前單純依賴衛(wèi)星數(shù)據(jù)構(gòu)建的全球巖石圈磁場(chǎng)球諧模型可靠階數(shù)可達(dá)90階(Olsen et al., 2017)，對(duì)應(yīng)空間分辨率為2°(在赤道地區(qū)約223 km).因此，采用這些模型的可靠波段對(duì)區(qū)域航磁異常匯編數(shù)據(jù)的長(zhǎng)波成分進(jìn)行評(píng)價(jià)與替換是提高區(qū)域磁力異常數(shù)據(jù)匯編質(zhì)量的關(guān)鍵步驟之一.顯然，構(gòu)建“全波段”均可靠的區(qū)域磁力異常匯編數(shù)據(jù)，作為一個(gè)國(guó)家或地區(qū)的基礎(chǔ)地球物理場(chǎng)數(shù)據(jù)，將能更好地服務(wù)于教學(xué)、科研以及各個(gè)行業(yè)和領(lǐng)域的廣泛應(yīng)用(例如輔助導(dǎo)航, 鄭暉等, 2012).

對(duì)于全球數(shù)據(jù)匯編，可以首先對(duì)其進(jìn)行球諧展開，然后將相應(yīng)長(zhǎng)波波段的球諧系數(shù)采用基于衛(wèi)星磁測(cè)數(shù)據(jù)構(gòu)建的全球巖石圈磁場(chǎng)模型的球諧系數(shù)替代(Maus et al., 2009; Lesur et al., 2016; 杜勁松等, 2017).對(duì)于區(qū)域數(shù)據(jù)匯編，盡管可以采用(超)長(zhǎng)剖面高精度磁測(cè)數(shù)據(jù)進(jìn)行調(diào)平處理，或者選擇編圖區(qū)內(nèi)部分切割線與測(cè)線組成若干磁場(chǎng)調(diào)平框架從而調(diào)整各個(gè)測(cè)區(qū)的磁場(chǎng)水平(尹航等, 2015)，但是更為高效的方法還是在長(zhǎng)波波段采用基于衛(wèi)星磁測(cè)數(shù)據(jù)構(gòu)建的全球巖石圈磁場(chǎng)模型預(yù)測(cè)數(shù)據(jù)進(jìn)行替換.

為了準(zhǔn)確且唯一地提取區(qū)域航磁異常數(shù)據(jù)的長(zhǎng)波成分，需要一種基函數(shù)完備、正則、調(diào)和且正交的區(qū)域巖石圈磁場(chǎng)建模方法，而修正球冠諧分析方法(Thébault et al., 2004, 2006a,b)正是滿足這些條件且是目前較為成熟的一種區(qū)域地磁場(chǎng)建模方法.該方法在地球外部一個(gè)球面錐體之內(nèi)對(duì)地磁場(chǎng)進(jìn)行建模，相比于傳統(tǒng)球冠諧分析方法(Haines, 1985a,b)，修正球冠諧分析方法考慮到兩組邊界條件而非一組，這兩組邊界條件引入了外源Legendre函數(shù)序列與Mehler函數(shù)序列，前者改善了衛(wèi)星高度數(shù)據(jù)的擬合與預(yù)測(cè)效果，后者改善了數(shù)據(jù)空缺區(qū)的上、下延拓效果.相比其它方法，R-SCHA方法具有計(jì)算量小、計(jì)算精度高的明顯優(yōu)勢(shì)，不僅可以融合地面、船載、航空與衛(wèi)星等多來源的地磁場(chǎng)不同分量(例如：總磁場(chǎng)強(qiáng)度、磁場(chǎng)三分量、磁場(chǎng)梯度等)的觀測(cè)數(shù)據(jù)，而且在分量轉(zhuǎn)換以及空間延拓等方面具有較高的計(jì)算精度(Ou et al., 2013; Vervelidou et al., 2018).因此，該方法已被廣泛應(yīng)用于區(qū)域性主磁場(chǎng)、巖石圈磁場(chǎng)、電離層與磁層磁場(chǎng)等的建模之中(徐文耀等, 2011; Talarn et al., 2017; Torta, 2020).

為了對(duì)航磁異常匯編數(shù)據(jù)中不可靠的低頻數(shù)據(jù)進(jìn)行提取、評(píng)價(jià)和替換，本文選用二維修正球冠諧分析(R-SCHA2D)方法(Thébault, 2008)，首先對(duì)航磁異常數(shù)據(jù)進(jìn)行R-SCHA2D反演建模，得到修正球冠諧系數(shù)，之后根據(jù)模型系數(shù)與波長(zhǎng)之間的關(guān)系提取長(zhǎng)波成分，最后使用全球巖石圈磁場(chǎng)模型的長(zhǎng)波成分對(duì)區(qū)域航磁異常匯編數(shù)據(jù)中的長(zhǎng)波成分進(jìn)行評(píng)價(jià)與替換.本文首先交代了R-SCHA2D的方法原理與計(jì)算算法，然后給出了區(qū)域航磁異常匯編數(shù)據(jù)長(zhǎng)波成分提取、評(píng)價(jià)與替換技術(shù)，最后通過仿真測(cè)試與實(shí)際應(yīng)用證明了該方法的有效性與實(shí)用性.

1 修正球冠諧分析方法

1.1 方法原理

修正球冠諧分析方法的建模區(qū)域Ω是由地心距為a和b的兩個(gè)曲面S1與S2以及以地心為頂點(diǎn)、半張角為θmax的球面椎體包圍而成(如圖1a所示).對(duì)于近似分布在同一高度上的數(shù)據(jù)(例如航磁異常匯編數(shù)據(jù))，建模區(qū)域上、下球冠面近似重合，此時(shí)修正球冠諧分析方法中Mehler函數(shù)的定義不再成立.針對(duì)此問題，Thébault(2008)提出了二維修正球冠諧分析方法，其建模范圍如圖1b所示.

當(dāng)建模區(qū)域的上、下球冠重合時(shí)，通過推導(dǎo)可以得到磁位異常在無源空間的級(jí)數(shù)表達(dá)形式(Thébault, 2008)：

(1)

圖1 三維(a)與二維(b)修正球冠諧分析的建模區(qū)域示意圖Fig.1 Sketch map of 3D (a) and 2D (b) modelling regions for revised spherical cap harmonic analysis

進(jìn)而，根據(jù)磁位異?？梢缘玫酱艌?chǎng)三分量異常：

(2)

航磁異常數(shù)據(jù)一般均為總磁場(chǎng)強(qiáng)度異常(ΔT)數(shù)據(jù)，而磁力異常矢量的模量一般遠(yuǎn)小于主磁場(chǎng)強(qiáng)度，因此可將總磁場(chǎng)強(qiáng)度異常視為磁場(chǎng)三分量異常在主磁場(chǎng)方向上的投影(Blakely, 1995; 管志寧, 2005; Hinze et al., 2013)，即

ΔT=ΔBxcosIcosD+ΔBycosIsinD+ΔBzsinI，(3)

其中，I與D分別為地球主磁場(chǎng)方向的傾角和偏角.

1.2 坐標(biāo)與分量轉(zhuǎn)換計(jì)算方法

在模型反演計(jì)算之前，首先需要將地理坐標(biāo)與地磁場(chǎng)三分量異常從球坐標(biāo)系轉(zhuǎn)換至球冠坐標(biāo)系(安振昌, 1992; De Santis et al., 1996)，坐標(biāo)與分量異常的轉(zhuǎn)換過程分為如下三步：

首先，通過式(4)將坐標(biāo)從球坐標(biāo)系(φ,θ,r)轉(zhuǎn)換到地心直角坐標(biāo)系(x,y,z)：

(4)

其次，通過式(5)對(duì)笛卡爾坐標(biāo)系進(jìn)行旋轉(zhuǎn)，使得笛卡爾坐標(biāo)系Z軸指向球冠中心：

(5)

其中，θ0與φ0為球冠中心的余緯度和經(jīng)度.最后，通過式(6)將坐標(biāo)從地心直角坐標(biāo)系轉(zhuǎn)換到球冠坐標(biāo)系中：

(6)

地磁場(chǎng)三分量異?？赏ㄟ^式(7)從球坐標(biāo)系轉(zhuǎn)換至球冠坐標(biāo)系：

(7)

其中

(8)

(9)

1.3 基于ΔT磁力異常數(shù)據(jù)的R-SCHA2D模型系數(shù)求解方法

式(1)與(2)的二維修正球冠諧級(jí)數(shù)表達(dá)式可被描述為一個(gè)線性系統(tǒng)：

(10)

其中，Gx、Gy、Gz分別為磁場(chǎng)三分量異常的二維修正球冠諧核函數(shù)矩陣，m為二維修正球冠諧系數(shù)矩陣.從式(3)可以看出，總磁場(chǎng)強(qiáng)度異?？梢暈榇艌?chǎng)三分量異常的線性組合，將式(7)代入式(10)中，首先將球冠坐標(biāo)系中地磁場(chǎng)三分量異常的核函數(shù)轉(zhuǎn)換至球坐標(biāo)系：

(11)

其次，根據(jù)式(3)對(duì)地磁場(chǎng)三分量異常的核函數(shù)進(jìn)行線性組合，即可得到基于總磁場(chǎng)強(qiáng)度異常數(shù)據(jù)的修正球冠諧系數(shù)解算方程：

=ΔT，(12)

上標(biāo)T表示矩陣的轉(zhuǎn)置.由于觀測(cè)數(shù)據(jù)不可避免地包含噪聲，因此將式(12)改寫為

(13)

其中數(shù)據(jù)加權(quán)矩陣為

(14)

式中εi為第i個(gè)ΔT磁力異常觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)差.式(13)為一個(gè)線性方程組，可以采用共軛梯度法進(jìn)行求解，但是在實(shí)際反演計(jì)算中，系數(shù)矩陣的條件數(shù)往往很大，會(huì)嚴(yán)重影響到共軛梯度算法的收斂性，因此本文采用預(yù)優(yōu)共軛梯度法(Pilkington, 1997)以改善系數(shù)矩陣的條件數(shù)，通過預(yù)處理使系數(shù)矩陣的特征值較集中分布，從而提高了迭代收斂速度及計(jì)算穩(wěn)定性.

1.4 航磁異常數(shù)據(jù)的補(bǔ)空與擴(kuò)邊方法

對(duì)于實(shí)際的航空磁力異常數(shù)據(jù)，往往存在數(shù)據(jù)空白區(qū)以及不規(guī)則的區(qū)域邊界，若不對(duì)其進(jìn)行處理而直接進(jìn)行建模，則在數(shù)據(jù)空白區(qū)可能構(gòu)建出虛假異常以及在邊界區(qū)域致使振蕩異常出現(xiàn).壓制這些現(xiàn)象最直接與最有效的方法即是對(duì)數(shù)據(jù)空區(qū)進(jìn)行補(bǔ)值以及對(duì)區(qū)域邊界外圍進(jìn)行擴(kuò)邊處理(安振昌, 2003).本文采用全球巖石圈磁場(chǎng)高階球諧模型，例如790階/次的EMM2017球諧模型(www.ngdc.noaa.gov/geomag/EMM/)或800階/次的WDMAMv2球諧模型(www.wdmam.org)，通過球諧解算得到與航磁異常匯編數(shù)據(jù)相同高度的ΔT磁力異常數(shù)據(jù)(陳康等, 2021)，從而完成對(duì)區(qū)域航磁異常匯編數(shù)據(jù)的空區(qū)補(bǔ)值和邊界區(qū)域擴(kuò)邊處理.

2 長(zhǎng)波成分的提取、評(píng)價(jià)與替換技術(shù)

對(duì)于區(qū)域航磁異常匯編數(shù)據(jù)長(zhǎng)波成分的提取、評(píng)價(jià)與替換，主要按照如下五個(gè)步驟進(jìn)行：

第一步，航磁異常數(shù)據(jù)預(yù)處理，主要包含三個(gè)方面：一是，對(duì)于實(shí)際航磁異常數(shù)據(jù)需要采用全球高階巖石圈磁場(chǎng)模型進(jìn)行空區(qū)補(bǔ)值與外圍擴(kuò)邊處理，在目標(biāo)區(qū)域外圍，若沒有足夠擴(kuò)邊的數(shù)據(jù)，則會(huì)由于邊界效應(yīng)將影響目標(biāo)區(qū)域內(nèi)部的建模精度；二是，實(shí)際航磁異常數(shù)據(jù)的分辨率可能較高，考慮到長(zhǎng)波成分提取的目標(biāo)任務(wù)以及計(jì)算量與計(jì)算效率問題，區(qū)域磁場(chǎng)模型無須擬合高頻成分，因此可以預(yù)先對(duì)補(bǔ)值與擴(kuò)邊之后的磁力異常數(shù)據(jù)進(jìn)行一定程度的降尺度和平滑處理，作為建模輸入磁力異常數(shù)據(jù)；三是，合理評(píng)估建模輸入磁力異常數(shù)據(jù)的噪聲，即每個(gè)測(cè)點(diǎn)磁力異常數(shù)據(jù)的標(biāo)準(zhǔn)差，若無法評(píng)價(jià)，則將數(shù)據(jù)加權(quán)矩陣每個(gè)對(duì)角線元素均設(shè)置為1 nT-1.

第三步，航磁異常長(zhǎng)波成分的提取與評(píng)價(jià).由于式(1)所示的基函數(shù)系具備正交性與完備性，因而修正球冠諧分析所得的各個(gè)系數(shù)是唯一的且可以確定的，而且非整數(shù)nk階與波長(zhǎng)l具有嚴(yán)格的對(duì)應(yīng)關(guān)系(Thébault et al., 2006b)，即

(15a)

(15b)

可見nk與球冠半張角θmax近似呈反比例關(guān)系，即球冠面積越小、nk越大，代表的波長(zhǎng)則越短.因此，通過截?cái)嗲蚬谥C級(jí)數(shù)的展開階數(shù)(如K階)而僅保留k≤K的系數(shù)，再通過正演計(jì)算即可完成長(zhǎng)波成分的提取.對(duì)于K值的確定，可以采用波長(zhǎng)與能量之間的譜分析通過對(duì)比航磁異常與衛(wèi)磁異常的能量譜曲線從而確定.但是為了簡(jiǎn)化計(jì)算，如引言所述，目前單純依賴衛(wèi)星數(shù)據(jù)構(gòu)建的全球巖石圈磁場(chǎng)球諧模型可靠階數(shù)可達(dá)90階(Olsen et al., 2017)，因此可以直接對(duì)比球諧90階及其以下的能量譜，即令l=90得到K值，從而判斷航磁異常長(zhǎng)波成分(即k≤K)的可靠性.修正球冠諧能量譜的計(jì)算采用Vervelidou和 Thébault (2015)提出的方法，其內(nèi)源與外源勒讓德函數(shù)序列在球冠面上的面積歸一化平均能量如式(16)所示：

(16)

第四步，航磁異常長(zhǎng)波成分的替換.通過上步對(duì)比航磁異常與衛(wèi)磁異常的能量譜曲線，確定截?cái)嗖ㄩL(zhǎng)，再采用衛(wèi)磁異常的長(zhǎng)波成分替換航磁異常的長(zhǎng)波成分.

第五步，航磁異常數(shù)據(jù)長(zhǎng)波成分替換之后的可靠性驗(yàn)證.由于衛(wèi)星測(cè)量高度較高(一般高于250 km)，距離巖石圈較遠(yuǎn)，衛(wèi)星磁力異常數(shù)據(jù)主要為長(zhǎng)波成分，因此可以采用替換了和未替換長(zhǎng)波成分的航磁異常匯編數(shù)據(jù)(ΔT)，分別與衛(wèi)星磁測(cè)數(shù)據(jù)聯(lián)合進(jìn)行修正球冠諧建模，進(jìn)而根據(jù)衛(wèi)星磁力異常數(shù)據(jù)的擬合差情況判斷航磁異常數(shù)據(jù)長(zhǎng)波成分替換的必要性，這是因?yàn)槿魞烧叩拈L(zhǎng)波成分不匹配則難以擬合觀測(cè)數(shù)據(jù).

3 仿真測(cè)試

為驗(yàn)證本文所提方法的有效性，本節(jié)首先使用全球巖石圈磁場(chǎng)模型通過正演計(jì)算仿真模擬航磁異常觀測(cè)數(shù)據(jù)與長(zhǎng)波成分理論數(shù)據(jù)，然后對(duì)模擬的航磁異常觀測(cè)數(shù)據(jù)進(jìn)行二維修正球冠諧建模進(jìn)而提取其長(zhǎng)波成分，最后將長(zhǎng)波成分的理論數(shù)據(jù)與提取數(shù)據(jù)進(jìn)行對(duì)比分析.由于澳大利亞航磁異常數(shù)據(jù)具有較高的分辨率與精度，因此在仿真測(cè)試與實(shí)際應(yīng)用中均采用澳大利亞作為研究區(qū)域.

3.1 仿真數(shù)據(jù)

考慮到計(jì)算量和計(jì)算效率，僅使用EMM2017全球巖石圈磁場(chǎng)模型(www.ngdc.noaa.gov/geomag/EMM/)的16到400階球諧系數(shù)，模擬計(jì)算了澳大利亞地區(qū)的總磁場(chǎng)強(qiáng)度異常數(shù)據(jù)作為測(cè)試的仿真觀測(cè)數(shù)據(jù)，如圖2a所示，數(shù)據(jù)分布范圍為東經(jīng)112°至東經(jīng)154°、南緯10°至南緯44°，網(wǎng)格間距為0.25°×0.25°，大地高度為500 m；同時(shí)，使用該模型的16～90階的正演計(jì)算數(shù)據(jù)作為長(zhǎng)波成分提取結(jié)果的檢核數(shù)據(jù)(圖2b).

圖2 仿真模擬的航磁異常觀測(cè)數(shù)據(jù)(a)與長(zhǎng)波成分理論數(shù)據(jù)(b)投影方式為等面積Hammer投影，中心經(jīng)線設(shè)置為東經(jīng)133°；色標(biāo)左端的NaN表示缺失數(shù)據(jù)區(qū)域.Fig.2 Simulated observation data (a) of aeromagnetic anomaly and theoretical data (b) of long wavelength componentsHammer equivalent projection is adopted. Central meridian is set to be 133°E. The NaN in the left side of the color bar denotes the data gap.

3.2 測(cè)試結(jié)果及其分析

將球冠的中心設(shè)置為東經(jīng)133°、南緯27°，球冠的半張角設(shè)置為27°；將Legendre函數(shù)序列的最大截?cái)嚯A數(shù)設(shè)為50階，Mehler函數(shù)序列的最大截?cái)嚯A數(shù)設(shè)為5階.對(duì)圖2a所示航磁異常觀測(cè)數(shù)據(jù)進(jìn)行二維修正球冠諧建模，所得模型預(yù)測(cè)數(shù)據(jù)與擬合殘差分別如圖3a與3b所示，提取出的小于與等于球諧90階的低頻成分?jǐn)?shù)據(jù)及其與檢核數(shù)據(jù)之間的殘差分別如圖3c與3d所示.

由圖3a與3b可以看出，模型較好地?cái)M合了觀測(cè)數(shù)據(jù)，僅少部分高頻異常數(shù)據(jù)受限于較低的展開階數(shù)而沒有得到擬合，一方面展開階數(shù)太高將大幅增加計(jì)算量，另一方面部分高頻信號(hào)的未擬合不會(huì)影響長(zhǎng)波成分的提取質(zhì)量.由圖3c與3d可以看出，從仿真數(shù)據(jù)中提取出的低頻成分在分布形態(tài)上與檢核數(shù)據(jù)一致，殘差的最大值為77.72 nT、最小值為-47.68 nT、平均值為0.26 nT、均方差為±7.74 nT，低頻成分的能量在分離中僅具有輕微的泄露，這說明本文提出的長(zhǎng)波成分提取技術(shù)具有有效性.

圖3 模型預(yù)測(cè)數(shù)據(jù)(a)及其擬合差(b)、提取的長(zhǎng)波成分?jǐn)?shù)據(jù)(c)及其誤差(d)Fig.3 Predicted data (a) and misfit (b) of the model, extracted data of long wavelength components (c) and its errors (d)

圖4為模擬航磁異常(EMM2017模型16～400階/次)與理論長(zhǎng)波成分(EMM2017模型16～90階/次)的球冠諧頻譜曲線.由于建模區(qū)域的半張角為27°，因此能譜曲線球冠諧階數(shù)的區(qū)間大小為180/27、初始階數(shù)為90/27.由圖4可以看出，理論長(zhǎng)波成分與航磁異常低于球諧90階的長(zhǎng)波成分的頻譜非常一致，這進(jìn)一步說明本文的長(zhǎng)波長(zhǎng)成分提取方法是可靠的，此外還說明對(duì)航磁異常數(shù)據(jù)進(jìn)行降尺度處理不會(huì)影響其長(zhǎng)波成分的提取.

圖4 模擬航磁異常與理論長(zhǎng)波成分的頻譜曲線Fig.4 Spectrum curves of simulated aeromagnetic and theoretical long-wavelength magnetic anomaly fields

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

為了進(jìn)一步說明本文所提方法在實(shí)際應(yīng)用中的有效性以及必要的數(shù)據(jù)處理步驟，本節(jié)將長(zhǎng)波長(zhǎng)成分提取、評(píng)價(jià)與替換技術(shù)應(yīng)用于澳大利亞地區(qū)的航磁異常匯編數(shù)據(jù).

4.1 澳大利亞航磁異常匯編數(shù)據(jù)

如圖5所示，本文所使用的航磁異常數(shù)據(jù)為第七代澳大利亞總磁場(chǎng)強(qiáng)度異常網(wǎng)格數(shù)據(jù)(Djomani et al., 2019).該網(wǎng)格數(shù)據(jù)的網(wǎng)格間距約為3弧秒(約80 m)，相較于第六代澳大利亞總磁場(chǎng)強(qiáng)度異常網(wǎng)格數(shù)據(jù)，新加入了234組航磁測(cè)量數(shù)據(jù)，且該網(wǎng)格數(shù)據(jù)已被調(diào)平到地形之上80 m的高度，網(wǎng)格數(shù)據(jù)所在的坐標(biāo)系為GDA94坐標(biāo)系.

圖5 澳大利亞航磁異常數(shù)據(jù)(a)與陸地高程(b)Fig.5 Aeromagnetic anomaly data (a) and land elevation (b) of theAustralia

4.2 數(shù)據(jù)處理

由于原始數(shù)據(jù)的數(shù)據(jù)量較大，在局部區(qū)域具有上萬(wàn)nT的異常值存在，且存在數(shù)據(jù)空缺區(qū)(圖5a)，因此在使用修正球冠諧分析方法建模之前，需要對(duì)原始數(shù)據(jù)進(jìn)行以下四步的處理：

(1)考慮到計(jì)算量，首先將原始數(shù)據(jù)降尺度為0.05°×0.05°間隔的網(wǎng)格數(shù)據(jù)；

(2)使用EMM2017全球巖石圈磁場(chǎng)模型的16～790階，對(duì)數(shù)據(jù)空區(qū)進(jìn)行補(bǔ)值，如圖6a所示；

(3)經(jīng)過降采樣處理后，部分高幅值異常變成了奇異點(diǎn)，故對(duì)降采樣之后的數(shù)據(jù)再進(jìn)行基于滑動(dòng)球冠的高斯濾波，高斯濾波器的濾波半徑設(shè)置為100 km，滑動(dòng)球冠的半張角為0.25°，對(duì)高斯濾波后的0.05°×0.05°網(wǎng)格間距數(shù)據(jù)(圖6b)進(jìn)行均值降采樣到0.25°×0.25°網(wǎng)格間距(圖6c)，并計(jì)算了各個(gè)網(wǎng)格點(diǎn)的標(biāo)準(zhǔn)差(圖6d)，作為反演建模的觀測(cè)數(shù)據(jù)誤差；

(4)將數(shù)據(jù)從GDA94大地坐標(biāo)系轉(zhuǎn)換到球冠坐標(biāo)系.

圖6 澳大利亞航磁異常數(shù)據(jù)處理結(jié)果(a) 降采樣為0.5°×0.5°網(wǎng)格間距的數(shù)據(jù)擴(kuò)邊結(jié)果； (b) 高斯平滑結(jié)果； (c) 高斯平滑濾除的高頻異常；(d) 0.25°×0.25°滑動(dòng)窗口統(tǒng)計(jì)計(jì)算得到的高頻異常的標(biāo)準(zhǔn)差分布.Fig.6 Processed results of aeromagnetic anomaly data of the Australia(a) Extended result based on the downscaled grid data with resolution of 0.5°×0.5°; (b) The result by Gaussian smooth; (c) The filtered high-frequency anomalies by Gaussian smooth; (d) Standard deviation map of the filtered high-frequency anomalies by statistically calculation in a moving window with a size of 0.25°×0.25°.

4.3 長(zhǎng)波成分的提取、評(píng)價(jià)與替換

對(duì)圖6b所示航磁異常數(shù)據(jù)進(jìn)行二維修正球冠諧分析，參數(shù)設(shè)置與仿真測(cè)試的相同，提取出的小于和等于球諧90階的長(zhǎng)波成分與短波成分?jǐn)?shù)據(jù)分別如圖7a與7b所示.前已述及，由于沒有觀測(cè)數(shù)據(jù)的約束，在區(qū)域磁場(chǎng)建模時(shí)，數(shù)據(jù)空區(qū)會(huì)產(chǎn)生虛假異常，進(jìn)而在處理轉(zhuǎn)換中引入較大誤差，故對(duì)數(shù)據(jù)空區(qū)進(jìn)行補(bǔ)值是有必要的，但是由于不同數(shù)據(jù)集之間的差異，補(bǔ)值同樣也會(huì)在數(shù)據(jù)中引入錯(cuò)誤的低頻成分，這可以明顯地從圖6a中看出，在數(shù)據(jù)縫合線兩邊的數(shù)據(jù)存在明顯的系統(tǒng)差異，從分離出的長(zhǎng)波長(zhǎng)成分(圖7a)中也可以看出，對(duì)數(shù)據(jù)進(jìn)行補(bǔ)值也引入了錯(cuò)誤的長(zhǎng)波長(zhǎng)成分.除了在補(bǔ)值區(qū)域與航磁異常數(shù)據(jù)之間的縫合區(qū)域之外，對(duì)比圖2b與圖7a可知，由航磁數(shù)據(jù)提取的長(zhǎng)波長(zhǎng)成分與EMM2017全球巖石圈磁場(chǎng)模型中的低頻成分(16～90階)的空間分布形態(tài)基本一致，但是磁力異常的幅值變化存在較大差異.

圖7 提取的澳大利亞航磁異常的長(zhǎng)波成分(a)與短波成分(b)Fig.7 Extracted data sets of long (a) and short (b) wavelength components from aeromagnetic anomaly data of the Australia

由航磁數(shù)據(jù)提取的長(zhǎng)波長(zhǎng)成分(圖7a)與EMM2017全球巖石圈磁場(chǎng)模型中的低頻成分(16～90階)(圖2b)在幅值方面存在較大差異.但是，該差異主要體現(xiàn)在航磁異常數(shù)據(jù)區(qū)域，而在補(bǔ)值與擴(kuò)邊區(qū)域差異微弱，這從側(cè)面說明計(jì)算方法是正確的、計(jì)算結(jié)果是可靠的.EMM2017全球巖石圈磁場(chǎng)模型是在EMAG2v3全球網(wǎng)格數(shù)據(jù)(Meyer et al., 2017)的基礎(chǔ)上進(jìn)行球諧分析構(gòu)建的，其16～133階采用了MF7模型(MF6模型的升級(jí)版, Maus et al., 2008)進(jìn)行了替換.為了說明MF7模型16～90階的可靠性，本文又采用最新的全球巖石圈磁場(chǎng)模型GRIMM_L120(Lesur et al., 2013)的16～90階、LCS-1(Olsen et al., 2017)與CM6(Sabaka et al., 2020)的16～90階、CHAOS7.6(Finlay et al., 2020)的21～90階球諧系數(shù)解算了澳大利亞地區(qū)航空高度的磁力異常，如圖8所示，其統(tǒng)計(jì)參數(shù)見表1.通過比較可以看出，五個(gè)模型解算的長(zhǎng)波磁力異常無論是在空間分布形態(tài)還是幅值變化方面均不存在較大的差異，這說明圖7a所示的航磁異常長(zhǎng)波成分自身是不可靠的.

圖8 不同全球巖石圈磁場(chǎng)模型解算的澳大利亞磁力異常長(zhǎng)波成分(a) GRIMM_L120模型的16～90階； (b) LCS-1模型的16～90階； (c) CM6模型的16～90階； (d) CHAOS7.6模型的21～90階.Fig.8 Magnetic anomaly data sets with long wavelength components over the Australia from different global lithospheric magnetic field modelsGRIMM_L120 (a), LCS-1 (b) and CM6 (c) models with spherical harmonic (SH) degrees 16～90,and CHAOS7.6 (d) model with SH degrees 21～90.

圖9為衛(wèi)磁異常、原始航磁異常與替換過低頻的航磁異常的頻譜曲線，可見未替換過低頻的航磁異常與衛(wèi)磁異常在球冠諧30階以下的低頻成分存在較大差異，鑒于衛(wèi)磁異常的可靠階數(shù)可達(dá)90階，因此證明原始航磁異常中的低頻成分是不可靠的，尤其是30階以下的波長(zhǎng)成分.在90階以上，單純基于衛(wèi)星磁測(cè)數(shù)據(jù)構(gòu)建的全球巖石圈磁場(chǎng)模型(如GRIMM_L120、LCS-1、CM6和CHAOS7.6)由于空間分辨率有限導(dǎo)致其能量低于EMM2017和航磁異常的能量，EMM2017模型由于聯(lián)合了衛(wèi)磁和航磁異常數(shù)據(jù)因而其能譜曲線與航磁異常的能譜曲線比較一致，而在200階以上，EMM2017的能量略高于航磁異常的能量，這是由于在提取航磁異常長(zhǎng)波成分時(shí)對(duì)原始數(shù)據(jù)進(jìn)行了降尺度和高斯平滑處理，但是這并不影響航磁異常長(zhǎng)波成分提取的可靠性.此外，通過圖9還可以發(fā)現(xiàn)，在所有單純基于衛(wèi)星磁測(cè)數(shù)據(jù)構(gòu)建的全球巖石圈磁場(chǎng)模型中，CHAOS7.6與LCS-1模型與實(shí)際航磁異常數(shù)據(jù)的一致性最好，尤其是在球冠諧110階以上的波長(zhǎng)成分，但是CHAOS7.6模型未給出16～20階的巖石圈磁場(chǎng)球諧系數(shù)，因而其能量在球冠諧30階以下偏低.

圖9 航磁異常與衛(wèi)磁異常的頻譜曲線Fig.9 Spectrum curves of aeromagnetic and satellite magnetic anomaly fields

此外，由表1可以看出，航磁異常提取的長(zhǎng)波成分的平均值與全球巖石圈磁場(chǎng)模型的平均值差異較大，這是由于這些差異主要來源于長(zhǎng)波成分，即球冠諧30階以下的波長(zhǎng)部分.

由于CM6模型(Sabaka et al., 2020)采用了?rsted、SAC-C、CHAMP與Swarm大量的衛(wèi)星磁測(cè)數(shù)據(jù)，并且采用了綜合法進(jìn)行地磁場(chǎng)建模，各種起源的地磁場(chǎng)之間的耦合性更強(qiáng)，其主磁場(chǎng)與變化磁場(chǎng)也被廣泛用于磁力異常數(shù)據(jù)的歸算處理，因此綜合考慮之后選取CM6模型的16～90階球諧系數(shù)解算的磁力異常數(shù)據(jù)對(duì)提取的澳大利亞航磁異常長(zhǎng)波數(shù)據(jù)進(jìn)行替換.具體計(jì)算步驟為：

(1)采用構(gòu)建的二維修正球冠諧模型解算澳大利亞原始航磁異常數(shù)據(jù)網(wǎng)格節(jié)點(diǎn)上的長(zhǎng)波成分磁力異常，將其從原始航磁異常數(shù)據(jù)之中扣除，得到剩余的短波成分磁力異常；

表1 澳大利亞總磁場(chǎng)強(qiáng)度異常長(zhǎng)波成分統(tǒng)計(jì)參數(shù)表Table 1 Statistic parameters of magnetic anomaly data sets with long wavelength components over the Australia

(2)采用CM6模型的16～90階球諧系數(shù)解算澳大利亞原始航磁異常數(shù)據(jù)網(wǎng)格節(jié)點(diǎn)上的長(zhǎng)波成分磁力異常，將其加入步驟(1)長(zhǎng)波成分提取之后的剩余短波成分磁力異常，即可得到最終的航磁異常數(shù)據(jù)(圖10a)；

(3)若需要對(duì)航磁異?？諈^(qū)進(jìn)行補(bǔ)值以及對(duì)外圍進(jìn)行擴(kuò)邊，則可以將EMM2017模型的16～90階球諧系數(shù)采用CM6模型的16～90階球諧系數(shù)進(jìn)行替換，得到新的16～790階的全球巖石圈磁場(chǎng)球諧系數(shù)之后，再根據(jù)空區(qū)待補(bǔ)值與待擴(kuò)邊數(shù)據(jù)的坐標(biāo)進(jìn)行磁力異常解算，將其與第(2)步處理結(jié)果組合，即可得到空區(qū)補(bǔ)值與擴(kuò)邊之后的航磁異常格網(wǎng)數(shù)據(jù)(圖10b).

進(jìn)一步地，為了驗(yàn)證航磁異常匯編數(shù)據(jù)長(zhǎng)波成分的替換效果，分別采用替換了和未替換長(zhǎng)波成分的航磁異常匯編數(shù)據(jù)(ΔT)，與CHAMP衛(wèi)星磁測(cè)數(shù)據(jù)(ΔBx、ΔBy與ΔBz)聯(lián)合進(jìn)行修正球冠諧建模.其中，所采用的CHAMP衛(wèi)星磁測(cè)數(shù)據(jù)如圖11所示，其處理流程和相關(guān)參數(shù)與全球巖石圈磁場(chǎng)模型GRIMM_L120(Lesur et al., 2013)的一致.采用相同的建模參數(shù)，基于替換了長(zhǎng)波成分的航磁異常匯編數(shù)據(jù)與CHAMP衛(wèi)星磁測(cè)數(shù)據(jù)聯(lián)合建模預(yù)測(cè)的衛(wèi)星測(cè)點(diǎn)處磁力異常及其擬合差如圖12所示，而基于未替換長(zhǎng)波成分的航磁異常匯編數(shù)據(jù)與CHAMP衛(wèi)星磁測(cè)數(shù)據(jù)聯(lián)合建模預(yù)測(cè)的衛(wèi)星測(cè)點(diǎn)處磁力異常及其擬合差如圖13所示，兩種情況之下的相關(guān)統(tǒng)計(jì)參數(shù)見表2與表3.由此可以看出，替換了長(zhǎng)波成分的航磁異常匯編數(shù)據(jù)與衛(wèi)星磁測(cè)數(shù)據(jù)融合效果更好，而由航磁異常匯編數(shù)據(jù)提取的長(zhǎng)波成分與衛(wèi)星磁測(cè)數(shù)據(jù)存在不相容性或非一致性.

但是，從圖12和表3可以看出，即使替換了航磁異常長(zhǎng)波成分，在衛(wèi)星磁測(cè)數(shù)據(jù)擬合殘差中還存在一些趨勢(shì)信號(hào)，筆者多次調(diào)試航磁異常與衛(wèi)磁異常的擬合權(quán)重，但是依然無法擬合這些殘存信號(hào).Vervelidou等(2018)認(rèn)為這些殘余信號(hào)來源于衛(wèi)磁數(shù)據(jù)處理中未完全消除的外源場(chǎng)(例如電離層和磁層)信號(hào)，且分量之間在物理上存在不耦合性，因此在建模時(shí)無法完全擬合.為了驗(yàn)證這種觀點(diǎn)，筆者又計(jì)算了GRIMM_L120、LCS-1、CM6和CHAOS7.6模型對(duì)衛(wèi)磁數(shù)據(jù)的擬合情況，發(fā)現(xiàn)這些殘余信號(hào)的確存在，如圖14所示.對(duì)比圖12與圖14可以發(fā)現(xiàn)，本文所構(gòu)建模型的擬合效果更好，這是由于航磁異常90階以上的中短波長(zhǎng)成分的可靠性更高.除此之外，筆者認(rèn)為這些無法擬合的趨勢(shì)信號(hào)還有可能來源于主磁場(chǎng)與巖石圈磁場(chǎng)之間的混疊效應(yīng).

表2 航磁數(shù)據(jù)與衛(wèi)星磁測(cè)數(shù)據(jù)聯(lián)合建模的模型預(yù)測(cè)數(shù)據(jù)統(tǒng)計(jì)參數(shù)表Table 2 Statistic parameters of model predictions by joint modelling of aeromagnetic data and satellite magnetic measurements

表3 航磁數(shù)據(jù)與衛(wèi)星磁測(cè)數(shù)據(jù)聯(lián)合建模的數(shù)據(jù)擬合殘差統(tǒng)計(jì)參數(shù)表Table 3 Statistic parameters of data misfits by joint modelling of aeromagnetic data and satellite magnetic measurements

5 結(jié)論

針對(duì)區(qū)域航磁異常匯編數(shù)據(jù)可能存在長(zhǎng)波成分的不可靠性，本文基于二維修正球冠諧分析(R-SCHA2D)以及單純基于衛(wèi)星磁測(cè)數(shù)據(jù)構(gòu)建的全球巖石圈磁場(chǎng)模型提出區(qū)域航磁異常匯編數(shù)據(jù)的長(zhǎng)波成分提取、評(píng)價(jià)與替換技術(shù).首先，使用EMM2017全球巖石圈磁場(chǎng)模型計(jì)算了仿真數(shù)據(jù)與檢核數(shù)據(jù)，仿真測(cè)試的結(jié)果證明了該方法可以正確地分離出數(shù)據(jù)中的低頻成分；然后，將本文提出的方法應(yīng)用于第七代澳大利亞總磁場(chǎng)強(qiáng)度異常網(wǎng)格數(shù)據(jù)，采用單純由衛(wèi)星磁測(cè)數(shù)據(jù)構(gòu)建的可靠長(zhǎng)波成分替換了由數(shù)據(jù)補(bǔ)值與擴(kuò)邊以及原始航磁異常數(shù)據(jù)中自帶的一些不可靠的長(zhǎng)波成分，得到了各個(gè)頻段均更可靠的澳大利亞及鄰區(qū)的航磁異常數(shù)據(jù)；最后，采用替換了和未替換長(zhǎng)波成分的航磁異常匯編數(shù)據(jù)，與CHAMP衛(wèi)星磁測(cè)三分量數(shù)據(jù)聯(lián)合進(jìn)行修正球冠諧建模，結(jié)果顯示，替換了長(zhǎng)波成分的航磁異常匯編數(shù)據(jù)與衛(wèi)星磁測(cè)數(shù)據(jù)融合效果更好，而由航磁異常匯編數(shù)據(jù)提取的長(zhǎng)波成分與衛(wèi)星磁測(cè)數(shù)據(jù)存在不相容性或非一致性.仿真測(cè)試與實(shí)際應(yīng)用均證明了本文所提的區(qū)域航磁異常數(shù)據(jù)長(zhǎng)波成分提取、評(píng)價(jià)與替換技術(shù)具有有效性和實(shí)用性，因此可以在區(qū)域航磁異常數(shù)據(jù)匯編中發(fā)揮重要作用.

圖10 最終合成的澳大利亞航磁異常數(shù)據(jù)(a) 空區(qū)未補(bǔ)值與未擴(kuò)邊； (b) 空區(qū)補(bǔ)值與擴(kuò)邊.Fig.10 Finally compiled aeromagnetic anomaly data of the Australia(a) The result of no data supplement and no extension; (b) The result after data supplement and extension.

圖13 使用未替換過低頻成分的航磁數(shù)據(jù)與衛(wèi)星數(shù)據(jù)聯(lián)合建模的模型預(yù)測(cè)數(shù)據(jù)與擬合殘差(a)與(b)分別為航磁異常(ΔT)預(yù)測(cè)數(shù)據(jù)與擬合殘差； (c)與(d)分別為衛(wèi)星磁測(cè)北向分量(ΔBx)的預(yù)測(cè)數(shù)據(jù)與擬合殘差； (e)與(f)分別為衛(wèi)星磁測(cè)東向分量(ΔBy)的預(yù)測(cè)數(shù)據(jù)與擬合殘差； (g)與(h)分別為衛(wèi)星磁測(cè)徑向分量(ΔBz)的預(yù)測(cè)數(shù)據(jù)與擬合殘差.Fig.13 Predictions and data misfits of the model by joint modelling of aeromagnetic anomaly data whose long wavelength components are not replaced and satellite magnetic measurements(a) and (b) are predicted data and data misfit of aeromagnetic anomaly, respectively; (c) & (d), (e) & (f), and (g) & (h) are predicted data and data misfit of northern component (ΔBx), eastern component (ΔBy) and radial component (ΔBz) of satellite magnetic anomaly, respectively.

圖14 CM6模型對(duì)衛(wèi)星磁測(cè)數(shù)據(jù)的擬合差分布(a) 北向分量(ΔBx)； (b) 東向分量(ΔBy)； (c) 徑向分量(ΔBz).Fig.14 Data misfit distributions of satellite magnetic data by using CM6 model(a) Northern component (ΔBx), eastern component (ΔBy) and radial component (ΔBz).

致謝感謝Geoscience Australia提供澳大利亞航空磁力異常匯編數(shù)據(jù)！球面投影圖件采用了Generic Mapping Tools(GMT)軟件進(jìn)行繪制，兩位匿名審稿專家也為本文提出了寶貴的修改建議，在此一并表示衷心感謝！