利用全球地磁場模型研究ΔT與Tap之間的差異性特征

陳 康 胡正旺 杜勁松,3

1 廣西壯族自治區(qū)第七地質(zhì)隊(duì)，廣西壯族自治區(qū)柳州市柳堡路3號(hào)，545100 2 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，武漢市魯磨路388號(hào)，430074 3 中國地質(zhì)大學(xué)(武漢)地質(zhì)過程與礦產(chǎn)資源國家重點(diǎn)實(shí)驗(yàn)室，武漢市魯磨路388號(hào)，430074

將標(biāo)量磁力儀測量的總磁場強(qiáng)度T的大小減去主磁場T0的大小，即可得到總磁場強(qiáng)度異常(ΔT)。在實(shí)際數(shù)據(jù)處理與正反演及定量解釋時(shí)，往往將ΔT近似為磁力異常矢量Ta在主磁場方向的投影分量Tap[1-3]，此時(shí)ΔT就和其他分量異常一樣，與場源的磁化強(qiáng)度大小呈簡單的線性關(guān)系，而與主磁場T0無直接關(guān)系，并且在場源外部空間滿足拉普拉斯方程，從而具備調(diào)和性質(zhì)，進(jìn)而可以簡化其數(shù)據(jù)處理、正反演與定量解釋過程。但當(dāng)Ta幅值較大時(shí)，這種近似將會(huì)產(chǎn)生較大誤差。針對此誤差，部分學(xué)者[4-7]分別提出利用實(shí)測ΔT數(shù)據(jù)計(jì)算投影分量Tap的最優(yōu)化方法、等效源方法、迭代方法，并且其理論模型的實(shí)驗(yàn)和實(shí)際應(yīng)用均表明，在高磁環(huán)境下考慮將ΔT近似為Tap的誤差具有必要性。

但上述關(guān)于ΔT與Tap之間差異性問題的相關(guān)研究均是針對局部磁力異常在笛卡爾直角坐標(biāo)系中完成，而對于區(qū)域乃至全球性中、長波長ΔT與Tap磁力異常之間的差異程度、分布規(guī)律及影響因素研究尚未見到相關(guān)報(bào)道。因此，本文基于全球地磁場球諧模型，從不同波長范圍與不同觀測高度探討將ΔT近似為Tap的誤差的全球分布特征，旨在為全球巖石圈磁場模型的應(yīng)用以及基于地面、航空與衛(wèi)星ΔT數(shù)據(jù)構(gòu)建區(qū)域乃至全球巖石圈磁力異常場模型提供參考資料。

1 ΔT與Tap差異性的理論分析

如圖1所示，總磁場強(qiáng)度異常(ΔT)為總磁場強(qiáng)度(T)的大小與主磁場強(qiáng)度(T0)的大小之差，Tap為磁力異常矢量(Ta)在主磁場方向的投影分量，即

圖1 ΔT、Ta與Tap之間的關(guān)系示意圖

ΔT=|T|-|T0|

(1)

Tap=|Ta|cosθ

(2)

式中，θ為T0與Ta之間的夾角，其計(jì)算公式為：

(3)

根據(jù)文獻(xiàn)[8]理論進(jìn)行推導(dǎo)，ΔT與Tap之間的差異E為:

(4)

由于ΔT≤|Ta|，因此誤差E≥0。分析式(4)可知，當(dāng)ΔT為0時(shí)，誤差E達(dá)到最大值，即

Emax=|Ta|2/(2|T0|)

(5)

此時(shí)，|T|與|T0|等長且與|Ta|形成等腰三角形，T0與Ta之間的夾角θ′為：

(6)

袁曉雨等[8]認(rèn)為，誤差E與測點(diǎn)和場源之間的相對位置密切相關(guān)，誤差比例最大的測點(diǎn)位于磁力異常矢量Ta接近垂直于主磁場T0的區(qū)域，隨著Ta值增大，誤差按照近似指數(shù)關(guān)系迅速增大，且誤差最大值的極值主要分布在中緯度地區(qū)，對應(yīng)的磁力異常矢量Ta與主磁場T0的夾角主要分布在90°～120°范圍內(nèi)。但上述關(guān)于誤差E的幅值大小與空間分布特征以及影響因素的部分結(jié)論是基于局部平面直角坐標(biāo)系的理論模型所獲得，可能與實(shí)際情況存在偏離。另外，對于不同的波長范圍與觀測高度，誤差E的幅值大小及其變化規(guī)律是構(gòu)建與應(yīng)用全球巖石圈磁場模型的重要因素。因此，需要從實(shí)際情況出發(fā)，采用全球巖石圈磁場模型研究誤差E的幅值大小與空間分布特征以及影響因素。

2 基于全球巖石圈磁場模型的ΔT與Tap及其差異計(jì)算方法

為分析實(shí)際情況下ΔT與Tap之間的差異性特征，本節(jié)給出基于全球巖石圈磁場模型計(jì)算相關(guān)物理量的方法。

主磁場與巖石圈磁場分別起因于地球外核磁流體運(yùn)動(dòng)與巖石圈磁性結(jié)構(gòu)，兩者共同組成地磁場的內(nèi)源場，在地球表面及其以上空間，可以將內(nèi)源場的磁位表示為球諧級(jí)數(shù)展開形式。根據(jù)展開式可以計(jì)算球面局部指北直角參考系(x軸指北、y軸指東、z軸指向地心)中的磁場三分量(X、Y與Z)，其無奇異性球諧計(jì)算表達(dá)式參見文獻(xiàn)[9]。為避免地球扁率的影響，所有解算點(diǎn)位坐標(biāo)與參考坐標(biāo)系分別為基于WGS84參考橢球的大地坐標(biāo)(大地緯度B、大地經(jīng)度L、大地高度H)與橢球面局部指北坐標(biāo)系，坐標(biāo)與矢量旋轉(zhuǎn)方法參見文獻(xiàn)[10]。

根據(jù)地磁場的球諧能量譜特征[11]，一般將地磁場球諧展開式中1～15階作為主磁場，而將16階及其以上部分作為巖石圈磁場。因此，將相應(yīng)階數(shù)范圍的球諧系數(shù)代入球諧表達(dá)式中，可計(jì)算出主磁場三分量(X0、Y0與Z0)、巖石圈磁場三分量(ΔX、ΔY與ΔZ)與總磁場三分量(X、Y與Z)，進(jìn)而可得到|T|、|T0|、ΔT、|Ta|、Tap：

(7a)

(7b)

(7c)

(7d)

(7e)

將式(7)代入式(4)與式(5)即可分別得到將ΔT近似為Tap的誤差E與最大誤差Emax，將式(7)代入式(3)即可得到T0與Ta之間的夾角θ為：

(8)

3 結(jié)果與分析

采用EMM2017地磁場模型計(jì)算與分析實(shí)際情況下ΔT與Tap之間的差異性特征。該模型由最新的EMAG2v3[12]全球總磁場強(qiáng)度異常2′×2′網(wǎng)格數(shù)據(jù)通過球諧展開獲得，最高展開階數(shù)達(dá)790階次，同時(shí)包含2000～2022年期間1～15階主磁場及其長期變化。因此，采用該模型1～15階、16～790階分別作為主磁場與巖石圈磁場。

3.1 ΔT與Tap之間差異的大小與分布特征

首先采用EMM2017(16～790階)計(jì)算得到參考橢球面上的總磁場強(qiáng)度異常(ΔT)、ΔT與Tap之差(圖2))。從圖中可以看出，ΔT與Tap之差均為非負(fù)值，說明ΔT≥Tap，這與理論分析結(jié)論一致。此外，該差異在全球大部分地區(qū)均低于2 nT，僅在部分區(qū)域存在較大值。

圖2 ΔT及其與Tap之間的差異分布

由表1可以看出，對于參考橢球面與5 km高度面上16～790階的巖石圈磁場，該差異的最大值分別可達(dá)248.5 nT與120.2 nT。由于目前地面與航空磁測的精度均遠(yuǎn)高于該差異，因此在圖2(b)所示差異性較大地區(qū)的實(shí)際地面與航空ΔT數(shù)據(jù)處理、正演與反演時(shí)不能忽略該差異。此外，隨著全球巖石圈磁場模型展開階數(shù)的增高與計(jì)算高度的降低，ΔT與Tap之差的差異性會(huì)更加顯著，這說明在構(gòu)建或應(yīng)用高階全球巖石圈磁場模型時(shí)需要考慮該差異，但對于衛(wèi)星高度的ΔT磁力異常數(shù)據(jù)則無需考慮。

表1 ΔT、Tap及兩者之間的差異

3.2 ΔT與Tap之間差異的影響因素

由式(4)可知，誤差E與|T0|、|Ta|和ΔT均存在關(guān)系。圖3(a)與3(b)分別為采用EMM2017解算參考橢球面上的磁力異常模量|Ta|(16～790階)與主磁場強(qiáng)度模量之倒數(shù)(1/|T0|)(1～15階)。對比圖2(b)與圖3可以明顯看出，誤差E主要受|Ta|大小影響，而受主磁場強(qiáng)度|T0|影響很弱，并且誤差E未表現(xiàn)出隨緯度而變化的分布特征。

圖4(a)為誤差E(≥2.0 nT)與|Ta|2之間的統(tǒng)計(jì)關(guān)系圖，從圖中可以看出，在整體趨勢上誤差E與|Ta|2近似滿足線性關(guān)系(藍(lán)色虛線)。圖4(b)為E(≥2.0 nT)與|Ta|2/2|T0|之間的統(tǒng)計(jì)關(guān)系圖，從圖中可以看出，采用式(5)計(jì)算的誤差E確實(shí)為誤差最大值Emax(紅色虛線)。對比圖4(a)與4(b)也可以看出，誤差E受主磁場強(qiáng)度T0影響很弱，而圖中偏離直線的點(diǎn)受ΔT大小影響。進(jìn)一步基于圖2和圖3中ΔT、|Ta|與誤差E之間空間分布的相關(guān)性可知，誤差E受|Ta|影響最大、ΔT次之、|T0|最弱。

圖3 |Ta|與1/|T0|分布

圖4 誤差E與|Ta|2、|Ta|2/2|T0|的統(tǒng)計(jì)關(guān)系

3.3 ΔT與Tap之間的差異估計(jì)與壓制方法

根據(jù)式(5)，采用EMM2017地磁場模型(16～790階)計(jì)算得到Emax分布，結(jié)果如圖5(a)所示。對比圖5(a)與圖2(b)可以看出，兩者在空間分布上具有較好的一致性，兩者差異見圖5(b)。從圖中可以看出，Emax幅度大于實(shí)際誤差，因此采用式(5)計(jì)算E會(huì)過高估計(jì)誤差，但計(jì)算得到的結(jié)果一方面可以作為誤差估計(jì)的上限，另一方面可以用于判斷誤差較大的區(qū)域。

圖5 Emax及其與實(shí)際誤差E之間的差異分布

從式(4)也可以看出，在將Tap近似為實(shí)際ΔT數(shù)據(jù)進(jìn)行磁力異常三分量與模量轉(zhuǎn)換處理時(shí)，為減小誤差E，可以引入Ta趨近于ΔT的約束條件。通過表1也可以看出，在統(tǒng)計(jì)意義上，|Ta|與ΔT的幅值比較接近，但這僅對磁力異常矢量的大小進(jìn)行約束，并未約束其矢量方向，因此后文將分析誤差E與磁力異常矢量Ta方向之間的關(guān)系。

根據(jù)式(8)，采用EMM2017地磁場模型計(jì)算得到θ分布，結(jié)果如圖6(a)所示。對比圖6(a)與圖2(b)可以看出，誤差E與夾角θ無相關(guān)性。根據(jù)式(6)，采用EMM2017地磁場模型計(jì)算得到θ′分布，結(jié)果如圖6(b)所示，從圖中可以看出，誤差E與夾角θ′具有明顯的對應(yīng)關(guān)系，即當(dāng)E較小時(shí)，θ′接近90°；而當(dāng)E較大時(shí)，θ′逐漸增大(本次計(jì)算中夾角θ′未超過95°)。由式(6)可知，θ′受控于磁力異常模量|Ta|與主磁場強(qiáng)度模量|T0|之比值，這說明在對總磁場強(qiáng)度異常數(shù)據(jù)ΔT進(jìn)行轉(zhuǎn)換處理或磁場建模時(shí)，為減小將Tap近似為ΔT的誤差，可以對磁力異常矢量的方向進(jìn)行約束，即要求磁力異常矢量與主磁場垂直的分量的幅度較小。這點(diǎn)從圖1也可以看出，即欲使誤差E最小，則要求磁力異常矢量Ta與主磁場T0方向平行。

圖6 T0與Ta之間的夾角θ分布以及當(dāng)ΔT為0時(shí)T0與Ta之間的夾角θ′分布

Lesur等[13]在利用全球總磁場強(qiáng)度異常匯編ΔT數(shù)據(jù)反演巖石圈磁場球諧模型系數(shù)時(shí)采用了上述磁力異常矢量的方向約束，并且認(rèn)為引入該約束是為壓制Backus效應(yīng)[14]。但根據(jù)上述分析，本文認(rèn)為該約束應(yīng)該是為減小將Tap近似為ΔT的誤差，而非壓制Backus效應(yīng)。實(shí)際上，所謂Backus效應(yīng)，即由標(biāo)量磁場值反演磁場矢量具有多解性；而Tap方向即主磁場方向，為磁力異常矢量在主磁場方向的投影。因此，Backus效應(yīng)和Tap與ΔT的差異性為兩個(gè)不同的問題。現(xiàn)今根據(jù)ΔT可以精確計(jì)算Tap[4-7]，因此在利用總磁場強(qiáng)度異常ΔT數(shù)據(jù)反演構(gòu)建巖石圈磁場模型時(shí)，無需引入磁力異常矢量的大小或方向約束。

4 結(jié) 語

本文采用全球地磁場模型EMM2017分別計(jì)算得到不同波長范圍和不同高度的ΔT與Tap及兩者差異E的全球分布，進(jìn)而分析該差異的幅值大小與空間分布特征以及影響差異性的因素，得到以下結(jié)論：

1)對不同波長范圍與不同高度情況下誤差E的幅值大小與空間分布特征進(jìn)行分析，可為構(gòu)建與應(yīng)用巖石圈磁場模型時(shí)是否需要考慮ΔT與Tap之間的差異性提供判斷依據(jù)。

2)在實(shí)際情況下，誤差E受|Ta|影響最大、ΔT次之、|T0|最弱，并且誤差E未表現(xiàn)出隨緯度而變化的分布特征，與T0和Ta之間的夾角也無相關(guān)性。

3)為壓制ΔT與Tap之間的差異性，可以引入|Ta|趨近于ΔT的約束條件，也可以將磁力異常矢量與主磁場垂直分量的幅度最小化作為約束，但目前研究表明，無論是等效源方法、最優(yōu)化方法還是迭代方法均能夠基于ΔT精確計(jì)算Tap，因此在利用總磁場強(qiáng)度異常ΔT數(shù)據(jù)反演構(gòu)建巖石圈磁場模型時(shí)，無需對磁力異常矢量的大小或方向進(jìn)行約束。

致謝：EMM2017全球地磁場模型來源于https:∥www.ngdc.noaa.gov/geomag/EMM/，在此表示衷心感謝！