區域航磁異常匯編數據的長波成分提取、評價與替換技術

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

1 中國地質大學(武漢)地球物理與空間信息學院， 武漢 430074 2 地球內部多尺度成像湖北省重點實驗室， 武漢 430074 3 地質過程與礦產資源國家重點實驗室， 武漢 430074

0 引言

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

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

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

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

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

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

1 修正球冠諧分析方法

1.1 方法原理

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

當建模區域的上、下球冠重合時，通過推導可以得到磁位異常在無源空間的級數表達形式(Thébault, 2008)：

(1)

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

進而，根據磁位異常可以得到磁場三分量異常：

(2)

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

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

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

1.2 坐標與分量轉換計算方法

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

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

(4)

其次，通過式(5)對笛卡爾坐標系進行旋轉，使得笛卡爾坐標系Z軸指向球冠中心：

(5)

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

(6)

地磁場三分量異常可通過式(7)從球坐標系轉換至球冠坐標系：

(7)

其中

(8)

(9)

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

式(1)與(2)的二維修正球冠諧級數表達式可被描述為一個線性系統：

(10)

其中，Gx、Gy、Gz分別為磁場三分量異常的二維修正球冠諧核函數矩陣，m為二維修正球冠諧系數矩陣.從式(3)可以看出，總磁場強度異常可視為磁場三分量異常的線性組合，將式(7)代入式(10)中，首先將球冠坐標系中地磁場三分量異常的核函數轉換至球坐標系：

(11)

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

=ΔT，(12)

上標T表示矩陣的轉置.由于觀測數據不可避免地包含噪聲，因此將式(12)改寫為

(13)

其中數據加權矩陣為

(14)

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

1.4 航磁異常數據的補空與擴邊方法

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

2 長波成分的提取、評價與替換技術

對于區域航磁異常匯編數據長波成分的提取、評價與替換，主要按照如下五個步驟進行：

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

第三步，航磁異常長波成分的提取與評價.由于式(1)所示的基函數系具備正交性與完備性，因而修正球冠諧分析所得的各個系數是唯一的且可以確定的，而且非整數nk階與波長l具有嚴格的對應關系(Thébault et al., 2006b)，即

(15a)

(15b)

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

(16)

第四步，航磁異常長波成分的替換.通過上步對比航磁異常與衛磁異常的能量譜曲線，確定截斷波長，再采用衛磁異常的長波成分替換航磁異常的長波成分.

第五步，航磁異常數據長波成分替換之后的可靠性驗證.由于衛星測量高度較高(一般高于250 km)，距離巖石圈較遠，衛星磁力異常數據主要為長波成分，因此可以采用替換了和未替換長波成分的航磁異常匯編數據(ΔT)，分別與衛星磁測數據聯合進行修正球冠諧建模，進而根據衛星磁力異常數據的擬合差情況判斷航磁異常數據長波成分替換的必要性，這是因為若兩者的長波成分不匹配則難以擬合觀測數據.

3 仿真測試

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

3.1 仿真數據

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

圖2 仿真模擬的航磁異常觀測數據(a)與長波成分理論數據(b)投影方式為等面積Hammer投影，中心經線設置為東經133°；色標左端的NaN表示缺失數據區域.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 測試結果及其分析

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

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

圖3 模型預測數據(a)及其擬合差(b)、提取的長波成分數據(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階/次)與理論長波成分(EMM2017模型16～90階/次)的球冠諧頻譜曲線.由于建模區域的半張角為27°，因此能譜曲線球冠諧階數的區間大小為180/27、初始階數為90/27.由圖4可以看出，理論長波成分與航磁異常低于球諧90階的長波成分的頻譜非常一致，這進一步說明本文的長波長成分提取方法是可靠的，此外還說明對航磁異常數據進行降尺度處理不會影響其長波成分的提取.

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

4 實際應用

為了進一步說明本文所提方法在實際應用中的有效性以及必要的數據處理步驟，本節將長波長成分提取、評價與替換技術應用于澳大利亞地區的航磁異常匯編數據.

4.1 澳大利亞航磁異常匯編數據

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

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

4.2 數據處理

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

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

(2)使用EMM2017全球巖石圈磁場模型的16～790階，對數據空區進行補值，如圖6a所示；

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

(4)將數據從GDA94大地坐標系轉換到球冠坐標系.

圖6 澳大利亞航磁異常數據處理結果(a) 降采樣為0.5°×0.5°網格間距的數據擴邊結果； (b) 高斯平滑結果； (c) 高斯平滑濾除的高頻異常；(d) 0.25°×0.25°滑動窗口統計計算得到的高頻異常的標準差分布.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 長波成分的提取、評價與替換

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

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

由航磁數據提取的長波長成分(圖7a)與EMM2017全球巖石圈磁場模型中的低頻成分(16～90階)(圖2b)在幅值方面存在較大差異.但是，該差異主要體現在航磁異常數據區域，而在補值與擴邊區域差異微弱，這從側面說明計算方法是正確的、計算結果是可靠的.EMM2017全球巖石圈磁場模型是在EMAG2v3全球網格數據(Meyer et al., 2017)的基礎上進行球諧分析構建的，其16～133階采用了MF7模型(MF6模型的升級版, Maus et al., 2008)進行了替換.為了說明MF7模型16～90階的可靠性，本文又采用最新的全球巖石圈磁場模型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階球諧系數解算了澳大利亞地區航空高度的磁力異常，如圖8所示，其統計參數見表1.通過比較可以看出，五個模型解算的長波磁力異常無論是在空間分布形態還是幅值變化方面均不存在較大的差異，這說明圖7a所示的航磁異常長波成分自身是不可靠的.

圖8 不同全球巖石圈磁場模型解算的澳大利亞磁力異常長波成分(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為衛磁異常、原始航磁異常與替換過低頻的航磁異常的頻譜曲線，可見未替換過低頻的航磁異常與衛磁異常在球冠諧30階以下的低頻成分存在較大差異，鑒于衛磁異常的可靠階數可達90階，因此證明原始航磁異常中的低頻成分是不可靠的，尤其是30階以下的波長成分.在90階以上，單純基于衛星磁測數據構建的全球巖石圈磁場模型(如GRIMM_L120、LCS-1、CM6和CHAOS7.6)由于空間分辨率有限導致其能量低于EMM2017和航磁異常的能量，EMM2017模型由于聯合了衛磁和航磁異常數據因而其能譜曲線與航磁異常的能譜曲線比較一致，而在200階以上，EMM2017的能量略高于航磁異常的能量，這是由于在提取航磁異常長波成分時對原始數據進行了降尺度和高斯平滑處理，但是這并不影響航磁異常長波成分提取的可靠性.此外，通過圖9還可以發現，在所有單純基于衛星磁測數據構建的全球巖石圈磁場模型中，CHAOS7.6與LCS-1模型與實際航磁異常數據的一致性最好，尤其是在球冠諧110階以上的波長成分，但是CHAOS7.6模型未給出16～20階的巖石圈磁場球諧系數，因而其能量在球冠諧30階以下偏低.

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

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

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

(1)采用構建的二維修正球冠諧模型解算澳大利亞原始航磁異常數據網格節點上的長波成分磁力異常，將其從原始航磁異常數據之中扣除，得到剩余的短波成分磁力異常；

表1 澳大利亞總磁場強度異常長波成分統計參數表Table 1 Statistic parameters of magnetic anomaly data sets with long wavelength components over the Australia

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

(3)若需要對航磁異常空區進行補值以及對外圍進行擴邊，則可以將EMM2017模型的16～90階球諧系數采用CM6模型的16～90階球諧系數進行替換，得到新的16～790階的全球巖石圈磁場球諧系數之后，再根據空區待補值與待擴邊數據的坐標進行磁力異常解算，將其與第(2)步處理結果組合，即可得到空區補值與擴邊之后的航磁異常格網數據(圖10b).

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

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

表2 航磁數據與衛星磁測數據聯合建模的模型預測數據統計參數表Table 2 Statistic parameters of model predictions by joint modelling of aeromagnetic data and satellite magnetic measurements

表3 航磁數據與衛星磁測數據聯合建模的數據擬合殘差統計參數表Table 3 Statistic parameters of data misfits by joint modelling of aeromagnetic data and satellite magnetic measurements

5 結論

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

圖10 最終合成的澳大利亞航磁異常數據(a) 空區未補值與未擴邊； (b) 空區補值與擴邊.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 使用未替換過低頻成分的航磁數據與衛星數據聯合建模的模型預測數據與擬合殘差(a)與(b)分別為航磁異常(ΔT)預測數據與擬合殘差； (c)與(d)分別為衛星磁測北向分量(ΔBx)的預測數據與擬合殘差； (e)與(f)分別為衛星磁測東向分量(ΔBy)的預測數據與擬合殘差； (g)與(h)分別為衛星磁測徑向分量(ΔBz)的預測數據與擬合殘差.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模型對衛星磁測數據的擬合差分布(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提供澳大利亞航空磁力異常匯編數據！球面投影圖件采用了Generic Mapping Tools(GMT)軟件進行繪制，兩位匿名審稿專家也為本文提出了寶貴的修改建議，在此一并表示衷心感謝！