王正濤， 張勇剛

（武漢大學(xué) 測(cè)繪學(xué)院，湖北 武漢 430079）

直到20年前，雖然人們對(duì)地球重力場(chǎng)的“靜態(tài)”部分了解相對(duì)準(zhǔn)確，但對(duì)其時(shí)間變化的了解僅限于較低的球諧階次。2002年，隨著GRACE（Gravity Recovery and Climate Experiment）任務(wù)的啟動(dòng)，情況得到了大幅改善；2018年，其后續(xù)任務(wù)GRACE-FO（Gravity Recovery and Climate Experiment-Follow On）成功實(shí)施，重力衛(wèi)星任務(wù)使人們能夠以前所未有的時(shí)空分辨率和精度測(cè)量地球重力隨時(shí)間的變化。

在年際和年代際時(shí)間尺度上，陸地水儲(chǔ)量和海洋質(zhì)量的變化似乎主導(dǎo)了觀測(cè)到的行星尺度重力變化。然而對(duì)于季節(jié)項(xiàng)的周期變化，觀測(cè)值和預(yù)測(cè)值擬合得并不理想。其原因是這些地表過(guò)程的模型仍然不完善，尤其是與陸地蓄水量有關(guān)的模型。同時(shí)，地球液核中的核心流在年際到10年的時(shí)間尺度上也能影響地球重力場(chǎng)變化，部分無(wú)法解釋的信號(hào)也可能與地核流動(dòng)有關(guān)［1］。地球流體核心的動(dòng)力學(xué)過(guò)程導(dǎo)致的質(zhì)量再分配主要以下述幾種不同的方式引起重力變化，如密度異常導(dǎo)致流體外核表層對(duì)流、分層密度層中的垂直運(yùn)動(dòng)產(chǎn)生密度遷移、由隨時(shí)間變化的核心流引起的壓力變化推動(dòng)和拉動(dòng)核幔邊界（core-mantle boundary， CMB）［2-3］。

最近的一些研究結(jié)果表明，觀測(cè)到的特定大地測(cè)量信號(hào)可能與地核動(dòng)力學(xué)過(guò)程有關(guān)。首先，Mandea等［4］證實(shí)，在對(duì)已知地表貢獻(xiàn)進(jìn)行校正后，衛(wèi)星記錄的時(shí)變重力場(chǎng)主要時(shí)空模態(tài)與磁場(chǎng)垂直分量的二階主要時(shí)空模態(tài)密切相關(guān)。研究涵蓋了2002年年中至2010年年中之間的有限時(shí)間跨度，在幾個(gè)地理位置，特別是在太平洋區(qū)域，剩余重力信號(hào)與磁場(chǎng)徑向分量的長(zhǎng)期加速度（徑向分量的二階導(dǎo)數(shù)）之間存在良好的相關(guān)性。此外，還對(duì)全球重力場(chǎng)和全球磁場(chǎng)徑向長(zhǎng)期加速度時(shí)間序列進(jìn)行了奇異值分解（singular value decomposition， SVD）。雖然每個(gè)時(shí)間序列的主導(dǎo)模式具有不同的空間結(jié)構(gòu)，但它們的年際變化具有良好的相關(guān)性［5］。同時(shí)后續(xù)研究［4-5］表明，良好的相關(guān)性一直持續(xù)到2015年底。文獻(xiàn)［6-7］研究探測(cè)到了全球定位系統(tǒng)（GPS）網(wǎng)絡(luò)記錄的地表位移中行星級(jí)尺度約為1mm的振幅信號(hào)，每個(gè)周期約為6年［6-7］。重力觀測(cè)和極移觀測(cè)［8］中也觀察到周期為6年的大地測(cè)量信號(hào)。在文獻(xiàn)［6-9］中，研究者認(rèn)為，地核中的動(dòng)力學(xué)過(guò)程是觀測(cè)信號(hào)的來(lái)源，主要是因?yàn)楹诵牧鲾y帶的軸向角動(dòng)量與6年期間的日長(zhǎng)變化之間存在聯(lián)系，這說(shuō)明通過(guò)大地測(cè)量觀測(cè)可以捕捉地核內(nèi)的動(dòng)力學(xué)過(guò)程［10］。

本文基于反演得到的GRACE重力場(chǎng)數(shù)據(jù)，利用多源數(shù)據(jù)模型在扣除各種非地核信號(hào)后，得到包含有地核信號(hào)的重力殘項(xiàng)序列，同時(shí)根據(jù)CHAOS-7地磁場(chǎng)模型，構(gòu)建出地磁場(chǎng)徑向分量、分量一階導(dǎo)數(shù)、分量二階導(dǎo)數(shù)的全球模型，并使用統(tǒng)計(jì)分解方法對(duì)重力場(chǎng)信號(hào)與磁場(chǎng)信號(hào)進(jìn)行處理，得到重力信號(hào)與磁場(chǎng)信號(hào)的空間模態(tài)與時(shí)變模態(tài)，發(fā)現(xiàn)重力信號(hào)第3模式的時(shí)變特征與磁場(chǎng)信號(hào)第2模式的時(shí)變信號(hào)相關(guān)性較好，在2014年以后，兩者均存在線性快速變化的特征，同時(shí)，其空間模態(tài)反映了太平洋低緯度地區(qū)較為明顯東西向“偶極子”形式的空間特征，與地磁場(chǎng)二階導(dǎo)數(shù)在太平洋地區(qū)的空間變化特征相對(duì)應(yīng)。

1 實(shí)驗(yàn)數(shù)據(jù)

1.1 GRACE重力場(chǎng)數(shù)據(jù)

GRACE重力場(chǎng)數(shù)據(jù)可以公開(kāi)獲得，主要由德克薩斯大學(xué)空間研究中心（CSR）、德國(guó)地學(xué)中心（GFZ）、美國(guó)噴氣推進(jìn)實(shí)驗(yàn)室（JPL）三大機(jī)構(gòu)發(fā)布。本文使用的是CSR GRACE/GRACE-FO RL06 Solutions數(shù)據(jù)，由CSR官網(wǎng)解算提供，其中，潮汐影響包括海洋潮汐、固體潮和極潮，都已在數(shù)據(jù)處理中去除，非潮汐的大氣和海洋的影響也同時(shí)在數(shù)據(jù)處理中剔除。選用CSR機(jī)構(gòu)提供的階數(shù)為60的模型，將官網(wǎng)發(fā)布的GRACE球諧系數(shù)解轉(zhuǎn)換為網(wǎng)格形式的質(zhì)量異常（詳細(xì)的計(jì)算方法及有關(guān)處理見(jiàn)GRACE CSR用戶手冊(cè)）。應(yīng)用的后處理包括：①添加對(duì)地球中心運(yùn)動(dòng)和地球扁率的獨(dú)立估計(jì)（C20改正）；②過(guò)濾以南北條紋為特征的相關(guān)誤差；③與地心有關(guān)的一階項(xiàng)改正；④采用ICE6G-D模型進(jìn)行冰川均衡調(diào)整（GIA）的更正。

1.2 GLDAS模型

全球陸地?cái)?shù)據(jù)同化系統(tǒng)（GLDAS）模型由美國(guó)宇航局協(xié)同美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心進(jìn)行發(fā)布。該模型使用的是衛(wèi)星觀測(cè)數(shù)據(jù)和地面觀測(cè)數(shù)據(jù)，通過(guò)先進(jìn)的陸地表面模型和數(shù)據(jù)同化技術(shù)可以得到地表流體的變化，并輸出陸地地表的各項(xiàng)物理參數(shù)。這一強(qiáng)大的新陸地表面建模系統(tǒng)集成了來(lái)自先進(jìn)觀測(cè)系統(tǒng)的數(shù)據(jù)，用以支持改進(jìn)水文氣象的研究與調(diào)查［11］。

根據(jù)GLDAS-NOAH模型，將總體陸地水（土壤濕度、根系水、植被水、徑流、冠層水、地表積雪等）儲(chǔ)量變化統(tǒng)一歸化成與GRACE數(shù)據(jù)相同的等效水高形式，并用到GRACE數(shù)據(jù)的觀測(cè)結(jié)果中，以扣除GRACE信號(hào)中陸地水儲(chǔ)量的成分，該模型時(shí)間分辨率為1個(gè)月，并以NetCDF4的格式按1°×1°的單元網(wǎng)格形式發(fā)布。

1.3 ECCO模型

海洋環(huán)流和氣候估算（ECCO）聯(lián)合會(huì)成立于1999年，主要對(duì)海洋環(huán)流及其氣候領(lǐng)域進(jìn)行詳盡的研究。如今，在NASA物理海洋學(xué)、建模構(gòu)建冰凍圈項(xiàng)目的持續(xù)支持下，ECCO正在提高海洋狀態(tài)的估計(jì)范圍和估計(jì)真實(shí)度。與通常的數(shù)據(jù)同化（如，在數(shù)值氣候預(yù)報(bào)中使用的數(shù)據(jù)同化）不同，ECCO提供了海洋隨時(shí)間演化的物理性特征一致性估計(jì)。ECCO的長(zhǎng)期目標(biāo)是綜合考慮海洋環(huán)流、海冰、生物化學(xué)、大氣等因素提供全球高分辨率的海洋狀態(tài)估計(jì)模型［12］。

使用ECCO Version 4、Release 4產(chǎn)品數(shù)據(jù)，ECCO V4R4綜合了總環(huán)流模型（MITgcm）和大多數(shù)可用的衛(wèi)星數(shù)據(jù)及現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)，以產(chǎn)生一個(gè)物理性一致的水循環(huán)閉合海洋模型，選取ECCO產(chǎn)品中海表面高度數(shù)據(jù)（SSH）以用于扣除海洋質(zhì)量變化信號(hào)。

1.4 WOA模型

世界海洋地圖集2018（WOA 2018）由美國(guó)國(guó)家海洋和大氣管理局氣候與全球變化計(jì)劃的撥款促成，該計(jì)劃由國(guó)家海洋數(shù)據(jù)中心（現(xiàn)為國(guó)家環(huán)境信息中心）統(tǒng)籌與管理，其目的是編制研究海洋學(xué)數(shù)據(jù)庫(kù)和基于這些數(shù)據(jù)庫(kù)進(jìn)行診斷研究。該地圖集所包含的數(shù)據(jù)收集于2018年世界海洋數(shù)據(jù)庫(kù)（WOD 2018）中，統(tǒng)計(jì)和客觀分析海洋的實(shí)地測(cè)量數(shù)據(jù)文件由國(guó)家環(huán)境信息中心（NCEI）在線分發(fā)。

實(shí)驗(yàn)采用的《2018年世界海洋地圖集》（WOA 2018）于2019年7月發(fā)布，其對(duì)早期版本的《世界海洋地圖集》進(jìn)行了數(shù)據(jù)更新，包括自上一次發(fā)布以來(lái)向世界海洋數(shù)據(jù)庫(kù)（WOD）添加的約300萬(wàn)個(gè)新海洋學(xué)模型以及更新了的質(zhì)量約束條件。發(fā)布的海洋物理參量有：溫度、鹽度、含氧量、硝酸鹽含量、磷酸鹽含量、硅酸鹽含量、電導(dǎo)率、密度。

1.5 ICE6G-D模型

目前國(guó)際上GIA全球模型有：ICE-3G、ICE-4G（VM2）、ICE-5G（VM2）以及ICE-6G_C（VM5a）等ICE模型和Paulson07模型和Geruo13模型；GIA區(qū)域模型主要針對(duì)南極地區(qū)，例如IJ05模型、IJ05_R2模型、W12和W12a模型。

多倫多大學(xué)根據(jù)ICE-NG（VMX）序列發(fā)布了冰川均衡調(diào)整模型，稱(chēng)為ICE-6G_C（VM5a），最初以球諧系數(shù)的形式發(fā)布，階數(shù)為256，并已被證明可為來(lái)自北美、歐亞大陸和南極洲的GPS觀測(cè)數(shù)據(jù)、GRACE衛(wèi)星反演得到的時(shí)變重力場(chǎng)數(shù)據(jù)以及放射性碳測(cè)定的相對(duì)海平面歷史數(shù)據(jù)提供精確的數(shù)據(jù)擬合，Purcell等［13］認(rèn)為ICE-6G_C模型存在改進(jìn)空間，提出了ICE-6G_D模型，在一定程度上優(yōu)化了模型結(jié)果，但仍然存在一定物理層面上的差異。

本文選用JPL發(fā)布的ICE-6G_D模型進(jìn)行GIA改正，以球諧系數(shù)的形式發(fā)布，階數(shù)為256階。

1.6 CHAOS-7模型

丹麥大學(xué)國(guó)家空間研究所基于低地球軌道衛(wèi)星Swarm、CryoSat-2、CHAMP、SAC-C和?rsted收集的磁場(chǎng)觀測(cè)以及地面觀測(cè)站月平均測(cè)量值的年度差異發(fā)布了1999年至2020年期間與時(shí)間相關(guān)的近地地磁場(chǎng)CHAOS-7模型。CHAOS-7模型包括一個(gè)球諧度為20的時(shí)變地核場(chǎng)、一個(gè)與LCS-1巖石圈場(chǎng)模型合并的25階以上的靜態(tài)內(nèi)部場(chǎng)、磁層場(chǎng)及其對(duì)應(yīng)感應(yīng)場(chǎng)的模型，描述衛(wèi)星矢量磁強(qiáng)計(jì)對(duì)準(zhǔn)的歐拉角估計(jì)數(shù)據(jù)以及CryoSat-2的磁強(qiáng)計(jì)校準(zhǔn)參數(shù)。只有來(lái)自磁靜區(qū)域的數(shù)據(jù)滿足嚴(yán)格的地磁平靜時(shí)間標(biāo)準(zhǔn)才用于磁場(chǎng)估計(jì)，使用的模型參數(shù)解算方法為迭代重加權(quán)正則化最小二乘法，與以前版本的CHAOS模型相比，在高階系數(shù)中，時(shí)間相關(guān)的地核內(nèi)源場(chǎng)的正則化條件更為寬松。

CHAOS-7是過(guò)去15年DTU開(kāi)發(fā)的一系列時(shí)間相關(guān)地磁場(chǎng)模型中的最新模型［14-16］。在過(guò)去的6年中，Swarm衛(wèi)星三星組提供的衛(wèi)星數(shù)據(jù)達(dá)到了全球完整和同質(zhì)的數(shù)據(jù)覆蓋。CHAOS-7的數(shù)據(jù)來(lái)源除了基于磁測(cè)衛(wèi)星上收集的磁強(qiáng)計(jì)數(shù)據(jù)，還有經(jīng)修訂的月平均地面觀測(cè)站長(zhǎng)期變化數(shù)據(jù)系列的更新版［15］。本次實(shí)驗(yàn)所提及的地磁場(chǎng)均為地核地磁場(chǎng)。

2 理論原理

2.1 數(shù)據(jù)處理流程

主要思路是利用多源數(shù)據(jù)模型對(duì)預(yù)處理完成后的GRACE重力信號(hào)進(jìn)行地表淺層物質(zhì)遷移的改正，以此來(lái)突出核幔邊界的重力信號(hào)。在對(duì)基于GRACE網(wǎng)格的重力時(shí)間序列進(jìn)行了所有影響（陸地蓄水、陸地冰和海洋的影響）的校正后，將PCAICA方法使用于剩余GRACE重力場(chǎng)以及地磁場(chǎng)模型，以研究其時(shí)空特異性。其主要的年際重力信號(hào)與磁場(chǎng)信號(hào)相關(guān)，含有發(fā)生在核幔邊界（CMB）的過(guò)程特征。流程示意圖如圖1。

圖1 實(shí)驗(yàn)流程示意Fig.1 Flowchart of experiment

2.2 GRACE時(shí)變重力場(chǎng)反演

基于不同的反演數(shù)據(jù)，主要有兩大類(lèi)反演方法：第一類(lèi)，基于GRACE的Level-2月地球重力場(chǎng)的球諧系數(shù)反演，這種方法首次由Wahr等［17］提出，是目前較為常用的方法，稱(chēng)為位理論法；第二類(lèi)，根據(jù)GRACE的Level-1B數(shù)據(jù)進(jìn)行反演，如質(zhì)量瘤法（mascon）法，該方法在一定程度上可以提高反演地表質(zhì)量變化的時(shí)空分辨率，但方法的技術(shù)原理較為復(fù)雜，處理流程難度較大。

在上述2種方法的基礎(chǔ)上，發(fā)展出新的方法，比如球諧系數(shù)擬合法、點(diǎn)質(zhì)量模型法和最優(yōu)先驗(yàn)?zāi)Ｐ头ǖ?。本文采用最為常用的位理論法?/p>

式中：Δh(θ，λ)為等效水高變化，反映地球地表水的質(zhì)量變化，其中θ為余緯，λ為經(jīng)度；a為地球平均參考半徑；ρa(bǔ)ve為地球平均密度；ρw為水的密度；Pˉlm(cosθ)為完全規(guī)格化的締合勒讓德函數(shù)，其中l(wèi)為階，m為次；kl為l階的負(fù)荷勒夫數(shù)；ΔCˉlm和ΔSˉlm為完全規(guī)格化的球諧系數(shù)變化，分別是cos (mλ)的系數(shù)、sin (mλ)的系數(shù)。式（1）即為用GRACE數(shù)據(jù)反演水儲(chǔ)量變化的基本公式。

2.3 地磁場(chǎng)模型數(shù)據(jù)

CHAOS-7模型的基本參數(shù)化與CHAOS系列早期版本的參數(shù)化相同，假設(shè)測(cè)量是在沒(méi)有電流的區(qū)域中收集的，因此，在電磁學(xué)的準(zhǔn)靜態(tài)近似下，矢量磁場(chǎng)B可以用標(biāo)量勢(shì)表示，使得B=-?V。磁標(biāo)量勢(shì)V=Vint+Vext，由內(nèi)部（主要是核心和巖石圈）源的勢(shì)Vint、外部（此處假設(shè)為磁層）源的勢(shì)Vext組成。內(nèi)部和外部都以球面諧波（SH）形式展開(kāi)。

對(duì)于內(nèi)源場(chǎng)，在地心地固系下（ECEF），對(duì)其進(jìn)行球諧展開(kāi)［16］，有

式中：(R，θ，λ)為球面坐標(biāo)；Nint為描述內(nèi)源場(chǎng)的最大階數(shù)；glm與hlm為描述內(nèi)源場(chǎng)的高斯系數(shù)，在CHOAS-7模型中，1～20階的內(nèi)源場(chǎng)高斯系數(shù)glm、hlm與時(shí)間有關(guān)，21階及以上的高斯系數(shù)是靜態(tài)的；Plm為施密特半歸一化締合勒讓德函數(shù)。

2.4 基于統(tǒng)計(jì)信息的信號(hào)分解

2.4.1 主成分分析（PCA）

將經(jīng)驗(yàn)正交函數(shù)分析方法（EOF）應(yīng)用于分解重力殘差項(xiàng)，用來(lái)研究改正后的GRACE重力場(chǎng)剩余信號(hào)和時(shí)空特性，該方法也稱(chēng)主成分分析（PCA），是一種分析矩陣數(shù)據(jù)中結(jié)構(gòu)特征、提取主要數(shù)據(jù)特征量的統(tǒng)計(jì)分解方法［18］。將一組GRACE單元格網(wǎng)點(diǎn)時(shí)間序列視為一組觀測(cè)值時(shí)間序列數(shù)據(jù)矩陣，對(duì)于整個(gè)格網(wǎng)的第i個(gè)格網(wǎng)點(diǎn)進(jìn)行m次觀測(cè)，則對(duì)應(yīng)該格網(wǎng)點(diǎn)的觀測(cè)向量Xi可表示為

由于實(shí)驗(yàn)地區(qū)有n個(gè)格網(wǎng)點(diǎn)，且對(duì)每個(gè)格網(wǎng)點(diǎn)都進(jìn)行過(guò)m次觀測(cè)，所以所有格網(wǎng)點(diǎn)上的距平觀測(cè)值時(shí)間序列按列排列可以用一個(gè)維度為s×n的矩陣Y表示，為

式中：s代表觀測(cè)時(shí)間序列的長(zhǎng)度。為了凸顯變化信息，式（4）對(duì)矩陣Y通常進(jìn)行距平處理，主成分分解的結(jié)果可以寫(xiě)為

式中：P為時(shí)變模態(tài)矩陣；E為空間模態(tài)矩陣。矩陣E中的s列數(shù)據(jù)代表觀測(cè)值Y的s個(gè)空間模態(tài)，即經(jīng)驗(yàn)正交函數(shù)，觀測(cè)值矩陣Y在第i個(gè)對(duì)應(yīng)空間模態(tài)上的投影pi=Yei就是該空間模態(tài)對(duì)應(yīng)的時(shí)變模態(tài)，也稱(chēng)為主成分，pi是矩陣P的第i列，ei是矩陣E的第i列。

式（5）右側(cè)的矩陣可以采用奇異值分解（SVD）方法進(jìn)行矩陣的統(tǒng)計(jì)信號(hào)分解，以研究在應(yīng)用所有校正后重力殘項(xiàng)仍然存在的信號(hào)及其時(shí)空變異性。

基于SVD分解可以將式（5）分解成為

式中：包含了正交化后的時(shí)變模態(tài)，即為包含按照從大到小順序排列的奇異值對(duì)角矩陣。

對(duì)于分解后保留的主要的r＜min (s，n)個(gè)空間模態(tài)和對(duì)應(yīng)的r個(gè)時(shí)變模態(tài)，觀測(cè)值矩陣可以通過(guò)式（7）進(jìn)行近似恢復(fù)：

式中：Yr為在保留前r個(gè)時(shí)變模態(tài)的情況下矩陣Y的近似矩陣；和Er的維度分別為s×r、r×r和n×r，且=Ir，ETrEr=Ir，每一個(gè)(j=1，2，3，…，r）被稱(chēng)為觀測(cè)值變量矩陣的“模式”，其中代表中的第j列數(shù)據(jù)，λj代表Λi對(duì)角線上第j個(gè)值，ej代表Er中第j列數(shù)據(jù)。

2.4.2 獨(dú)立主成分分析（ICA）

PCA方法是一種基于二階統(tǒng)計(jì)信息的分解方法，在分解過(guò)程中，不管是基于時(shí)間序列構(gòu)造的自協(xié)方差矩陣還是相關(guān)性矩陣，都只包含觀測(cè)信號(hào)中的二階統(tǒng)計(jì)信息，而忽略了觀測(cè)值概率密度函數(shù)的高階矩中所含的信號(hào)，且經(jīng)過(guò)PCA分解得到的不同模式僅表現(xiàn)為彼此不相關(guān)，而并非相互獨(dú)立。

基于上述討論，需要對(duì)傳統(tǒng)的PCA方法進(jìn)行改進(jìn)擴(kuò)展，為了從基于觀測(cè)數(shù)據(jù)的概率密度函數(shù)中獲取更多的信息，就需要在分解過(guò)程中顧及概率密度函數(shù)的高階統(tǒng)計(jì)信息作為影響因子，由此有關(guān)學(xué)者便提出獨(dú)立主成分分析（ICA）方法。ICA方法的基本思想是利用高階統(tǒng)計(jì)信息，使得分解得到的成分盡可能彼此獨(dú)立，旨在在沒(méi)有先驗(yàn)信息的條件下，利用高階統(tǒng)計(jì)信息，將觀測(cè)信號(hào)分解成多個(gè)彼此相互獨(dú)立的成分，可以視為是基于二階統(tǒng)計(jì)信息分解方法的一種拓展延伸。

將統(tǒng)計(jì)信息的信號(hào)分解問(wèn)題一般化，轉(zhuǎn)化為尋找將觀測(cè)值矩陣Y轉(zhuǎn)化為源信號(hào)矩陣S的函數(shù)f。

在多數(shù)應(yīng)用場(chǎng)景下，式（8）轉(zhuǎn)化關(guān)系為了簡(jiǎn)化計(jì)算，可以簡(jiǎn)單定義為線性轉(zhuǎn)換，因此，式（8）的函數(shù)變換可以用矩陣乘積的形式表示，為

式中：W為待求的解混矩陣。

目前常用的確定解混矩陣W的實(shí)用計(jì)算方法也是SVD分解方法，如式（10）：

ICA分解可以看作是2.4.1節(jié)中所講到的PCA分解的旋轉(zhuǎn)擴(kuò)展應(yīng)用，通過(guò)上述分析，基于正交模式旋轉(zhuǎn)的ICA分解可表示為

式中：Rr為維度是r×r的矩陣，為了得到ICA分解的模態(tài)，便需要確定此旋轉(zhuǎn)矩陣，該矩陣對(duì)Er和P^r進(jìn)行矩陣旋轉(zhuǎn)時(shí)使它們的列之間盡可能在統(tǒng)計(jì)上保持獨(dú)立，并且其旋轉(zhuǎn)矩陣本身也具有正交性，即RrRTr=Ir。那么，根據(jù)式（11）得到統(tǒng)計(jì)獨(dú)立源信號(hào)可表示為

基于式（12）得到的源信號(hào)的行在空間上盡可能達(dá)到獨(dú)立。以上分解稱(chēng)為ICA分解。

3 實(shí)驗(yàn)結(jié)果

3.1 重力場(chǎng)模型數(shù)據(jù)處理

為統(tǒng)一各個(gè)模型的時(shí)間跨度，分析2002年4月至2015年12月GRACE重力場(chǎng)的月時(shí)變序列，并對(duì)重力場(chǎng)模型進(jìn)行相應(yīng)的校正，提取與地核物理過(guò)程有關(guān)的深部地球重力信號(hào)，剔除的成分包括大氣負(fù)荷、陸地流域水文變化和海洋環(huán)流變化引起的海洋質(zhì)量變化以及GIA。

3.1.1 GRACE時(shí)變重力場(chǎng)

使用CSR機(jī)構(gòu)提供的GRACE重力場(chǎng)模型，該數(shù)據(jù)可以公開(kāi)獲得，并根據(jù)2.2節(jié)提及的相關(guān)原理，將其數(shù)據(jù)轉(zhuǎn)化為表征網(wǎng)格質(zhì)量異常的等效水高形式，由于GRACE的衛(wèi)星軌道形狀對(duì)重力場(chǎng)的低階項(xiàng)（C20項(xiàng)）不敏感，其低階帶諧項(xiàng)C20的精度相對(duì)較低，所以以衛(wèi)星激光測(cè)距技術(shù)（SLR）測(cè)得的C20項(xiàng)代替GRACE的C20項(xiàng)。得到的結(jié)果如圖2所示。

圖2 GRACE全球時(shí)變重力場(chǎng)等效水高（2002年4月）Fig.2 GRACE global time-varying gravity field equivalent water height(April 2022)

3.1.2 大氣負(fù)荷與海洋壓力

根據(jù)CSR發(fā)布的GRACE用戶手冊(cè)可知，對(duì)于全球大氣和海洋負(fù)載質(zhì)量變化所引起的重力場(chǎng)變化，在Level-1B階段，根據(jù)AOD1B（Atmosphere andOcean De-Aliasing Level-1B）Release 6數(shù)據(jù)進(jìn)行校正。

3.1.3 陸地水文成分

選用GLDAS水文模型來(lái)剔除GRACE重力場(chǎng)中的陸地水文成分，綜合總體陸地水儲(chǔ)量變化異常（土壤濕度、根系水、植被水、徑流、冠層水、地表積雪等），將其統(tǒng)一歸化成與GRACE數(shù)據(jù)相同的等效水高形式，調(diào)整空間分辨率并選取對(duì)應(yīng)的時(shí)間跨度用到GRACE數(shù)據(jù)的觀測(cè)結(jié)果中。

3.1.4 海洋質(zhì)量變化

實(shí)驗(yàn)使用ECCO Version 4、Release 4產(chǎn)品中發(fā)布的海洋表面高（SSH）數(shù)據(jù)，得到全球測(cè)高的總體海平面變化，聯(lián)合WOA18數(shù)據(jù)集的鹽度分層模型與溫度分層模型計(jì)算得到比容海平面變化，同時(shí)顧及冰后回彈對(duì)海面高的影響，計(jì)算得到由于全球海水質(zhì)量變化引起的等效水高［19］。

3.1.5 GIA影響

實(shí)驗(yàn)選用JPL發(fā)布的ICE-6G_D模型進(jìn)行GIA修正，將球諧系數(shù)統(tǒng)一轉(zhuǎn)化為等效水高的形式。

3.2 地磁場(chǎng)模型數(shù)據(jù)構(gòu)建

丹麥大學(xué)國(guó)家空間研究所使用CHAOS-7模型研究了地磁場(chǎng)的時(shí)間變化特征。本文基于CHAOS-7模型構(gòu)建地球表面以及核幔邊界處地磁場(chǎng)徑向分量強(qiáng)度和一階導(dǎo)數(shù)、二階導(dǎo)數(shù)模型。圖3—5中，F(xiàn)、F′、F″分別表示地磁場(chǎng)徑向強(qiáng)度、地磁場(chǎng)徑向強(qiáng)度關(guān)于時(shí)間的一階導(dǎo)數(shù)、地磁場(chǎng)徑向強(qiáng)度關(guān)于時(shí)間的二階導(dǎo)數(shù)。

圖3 全球表面地磁場(chǎng)徑向分量強(qiáng)度（2019年1月，截?cái)嚯A數(shù)20階）Fig.3 Radial magnetic field on surface of the earth(Jan.2019, up to SH degree 20)

圖4 全球表面地磁場(chǎng)徑向分量強(qiáng)度一階導(dǎo)數(shù)（2019年1月，截?cái)嚯A數(shù)20階）Fig.4 First time derivative of radial magnetic field on surface of the earth (Jan.2019, up to SH degree 20)

圖3—5顯示了時(shí)間歷元為2019年1月的CHAOS-7模型在地表的主磁場(chǎng)徑向分量、其一階導(dǎo)數(shù)（長(zhǎng)期變化，SV）和二階導(dǎo)數(shù)（持續(xù)加速度，SA）的全球分布圖，截?cái)嗟那蛑C階數(shù)為20階。

地球表面徑向磁場(chǎng)分量具有一些顯著特征，比如南大西洋的高緯度斑塊和赤道區(qū)域較弱的徑向磁場(chǎng)。徑向場(chǎng)一階導(dǎo)數(shù)顯示，最大的徑向場(chǎng)增加區(qū)域位于南美洲的東北角，同時(shí)有一條從南部非洲向西南延伸的徑向磁場(chǎng)帶。地球表面的徑向場(chǎng)二階導(dǎo)數(shù)顯示了太平洋中強(qiáng)烈的局部偶極結(jié)構(gòu)（在圖5的中部區(qū)域可見(jiàn)），東太平洋（包括夏威夷附近）出現(xiàn)負(fù)加速度、澳大利亞?wèn)|北部與西太平洋出現(xiàn)正加速度［20-21］。

圖5 全球表面地磁場(chǎng)徑向分量強(qiáng)度二階導(dǎo)數(shù)（2019年1月，截?cái)嚯A數(shù)20階）Fig.5 Second time derivative of radial magnetic field on surface of the earth (Jan.2019, up to SH degree 20)

圖6—8 顯示了時(shí)間歷元為2019 年1 月的CHAOS-7 模型在核幔邊界的主磁場(chǎng)徑向分量、其一階導(dǎo)數(shù)（長(zhǎng)期變化，SV）和二階導(dǎo)數(shù)（持續(xù)加速度，SA）的全球分布圖。針對(duì)核幔邊界的研究區(qū)域，相關(guān)研究者指出在幾年或更長(zhǎng)的時(shí)間尺度內(nèi)，地幔中電磁感應(yīng)在與時(shí)間有關(guān)的內(nèi)源場(chǎng)模型中起次要作用。同時(shí)，為了確保核幔邊界的地磁場(chǎng)強(qiáng)度（包括一階導(dǎo)數(shù)和二階導(dǎo)數(shù)）圖的圖像平滑與穩(wěn)定，實(shí)驗(yàn)選取的截?cái)嚯A數(shù)分別為13、17、15［16］。對(duì)于核幔邊界的信號(hào)而言，在非洲與南美洲之間的赤道地區(qū)，存在最強(qiáng)的地磁場(chǎng)徑向分量一階導(dǎo)數(shù)信號(hào)，該特征是區(qū)域徑向場(chǎng)向西運(yùn)動(dòng)的結(jié)果［15-16］；對(duì)于核幔邊界的二階導(dǎo)數(shù)圖，可以看出，位于東經(jīng)60°～90°的太平洋地帶存在大規(guī)模的加速現(xiàn)象，這與地表特征相對(duì)應(yīng)，但位于西經(jīng)60°～90°地區(qū)的大規(guī)模減速現(xiàn)象相對(duì)于地表特征來(lái)說(shuō)并不明顯。

圖6 核幔邊界地磁場(chǎng)徑向分量強(qiáng)度（2019年1月，截?cái)嚯A數(shù)13階）Fig.6 Radial magnetic field at core-mantle boundary (Jan.2019, up to SH degree 13)

圖7 核幔邊界地磁場(chǎng)徑向分量強(qiáng)度一階導(dǎo)數(shù)（2019年1月，截?cái)嚯A數(shù)17階）Fig.7 First time derivative of radial magnetic field at core-mantle boundary (Jan.2019, up to SH degree 17)

圖8 核幔邊界地磁場(chǎng)徑向分量強(qiáng)度二階導(dǎo)數(shù)（2019年1月，截?cái)嚯A數(shù)15階）Fig.8 Second time derivative of radial magnetic field at the core-mantle boundary (Jan.2019,up to SH degree 15)

3.3 重力場(chǎng)數(shù)據(jù)統(tǒng)計(jì)分解

以緯度-60°～60°之間的地區(qū)作為實(shí)驗(yàn)區(qū)域，采用PCA-ICA方法將實(shí)驗(yàn)區(qū)域內(nèi)的GRACE重力殘項(xiàng)序列分解為相應(yīng)的空間模態(tài)和時(shí)變模態(tài)。圖9展示了采用PCA-ICA分解后的不同數(shù)量GRACE模式的能量占比，圖中橫坐標(biāo)表示GRACE的模式數(shù)，縱坐標(biāo)表示模式能量之和占總能量的比例，例如橫坐標(biāo)為10時(shí)，其對(duì)應(yīng)點(diǎn)的縱坐標(biāo)表示前10個(gè)GRACE重力殘項(xiàng)模式的能量之和占所有GRACE模式能量之和的分?jǐn)?shù)。

圖9 不同數(shù)量GRACE模式的總能量占比Fig.9 Percentage of total energy in GRACE modes

根據(jù)圖9的統(tǒng)計(jì)結(jié)果，結(jié)合有關(guān)學(xué)者的研究成果，選取PCA-ICA分解的前6個(gè)模式作為主要模式進(jìn)行研究。圖10展示了全球?qū)嶒?yàn)區(qū)域6個(gè)主要模式的空間模態(tài)。對(duì)應(yīng)實(shí)驗(yàn)區(qū)域6個(gè)主要模式的時(shí)變模態(tài)見(jiàn)圖11所示。

圖10 PCA-ICA方法分解出的GRACE空間模態(tài)Fig.10 Spatial mode of GRACE using PCA-ICA

圖11 PCA-ICA方法分解出的前6個(gè)GRACE模式的時(shí)變模態(tài)Fig.11 Time-varying mode of GRACE modes using PCA-ICA

由10圖可知，前2個(gè)模式存在明顯的季節(jié)性特征，后續(xù)4個(gè)模式則反映了長(zhǎng)周期變化，其中除了第3模式，另外3個(gè)模式除了長(zhǎng)周期信號(hào)特征還反映了實(shí)驗(yàn)區(qū)域的季節(jié)性信號(hào)。前6個(gè)模式的能量占比分別為21.60%、13.23%、9.27%、5.45%、3.29%、2.72%。

3.4 地磁場(chǎng)數(shù)據(jù)統(tǒng)計(jì)分解

根據(jù)已有研究成果可知，重力殘項(xiàng)與磁場(chǎng)徑向分量的二階導(dǎo)數(shù)存在一定相關(guān)性，本實(shí)驗(yàn)重點(diǎn)針對(duì)地磁場(chǎng)模型徑向分量的二階導(dǎo)數(shù)數(shù)據(jù)，采用相同的PCA-ICA分解方法對(duì)實(shí)驗(yàn)區(qū)域（緯度-60°～60°）內(nèi)的磁場(chǎng)模型進(jìn)行統(tǒng)計(jì)分解。

3.4.1 地球表面地磁場(chǎng)

對(duì)于地磁場(chǎng)徑向分量的統(tǒng)計(jì)分解實(shí)驗(yàn)，與重力場(chǎng)數(shù)據(jù)相統(tǒng)一，選取PCA-ICA分解的前6個(gè)模式作為主要模式。不同數(shù)量地球表面地磁場(chǎng)模式的總能量占比見(jiàn)圖12。圖13展示全球?qū)嶒?yàn)區(qū)域6個(gè)主要模式的空間模態(tài)，前6個(gè)模式的信號(hào)能量占比分別為39.42%、17.73%、16.05%、8.53%、5.18%、4.09%。其對(duì)應(yīng)的6個(gè)模式的時(shí)變模態(tài)結(jié)果如圖14。

圖12 不同數(shù)量地球表面地磁場(chǎng)模式的總能量占比Fig.12 Percentage of total energy in geomagnetic field modes on surface of the earth

圖13 PCA-ICA方法分解出的地球表面地磁場(chǎng)（二階導(dǎo)數(shù)）空間模態(tài)Fig.13 Spatial mode of second time derivative of radial magnetic field on surface of the earth using PCA-ICA

圖14 PCA-ICA方法分解出的前6個(gè)地球表面地磁場(chǎng)（二階導(dǎo)數(shù)）時(shí)變模態(tài)Fig.14 Time-varying mode of second time derivative of radial magnetic field on surface of the earth using PCA-ICA

3.4.2 核幔邊界地磁場(chǎng)

對(duì)于核幔邊界地磁場(chǎng)徑向分量的統(tǒng)計(jì)分解實(shí)驗(yàn)，選取PCA-ICA分解的前6個(gè)模式作為主要模式。不同數(shù)量核幔邊界地磁場(chǎng)模式的總能量占比見(jiàn)圖15。圖16展示全球?qū)嶒?yàn)區(qū)域6個(gè)主要模式的空間模態(tài)，其信號(hào)能量占比分別為53.00%、17.70%、15.38%、6.72%、2.68%、2.03%。其對(duì)應(yīng)的6個(gè)模式的時(shí)變模態(tài)結(jié)果如圖17。

圖15 不同數(shù)量核幔邊界地磁場(chǎng)模式的總能量占比Fig.15 Percentage of total energy in geomagnetic field modes at the core-mantle boundary

圖16 PCA-ICA方法分解出的核幔邊界地磁場(chǎng)（二階導(dǎo)數(shù)）空間模態(tài)Fig.16 Spatial mode of second time derivative of radial magnetic field at the core-mantle boundary using PCA-ICA

圖17 PCA-ICA方法分解出的前6個(gè)地球表面地磁場(chǎng)（二階導(dǎo)數(shù)）模式時(shí)變模態(tài)Fig.17 Time-varying mode of second time derivative of radial magnetic field at the core-mantle boundary using PCA-ICA

3.5 重力信號(hào)與磁場(chǎng)信號(hào)的對(duì)比分析

根據(jù)圖11中6個(gè)模式對(duì)應(yīng)的時(shí)變模態(tài)信號(hào)特征分析可得，第1、4模式存在明顯的季節(jié)性周期特征，應(yīng)該是由于在進(jìn)行陸地水文信號(hào)、海洋質(zhì)量變化信號(hào)等非地核信號(hào)的扣除時(shí)并未扣除干凈，導(dǎo)致重力殘項(xiàng)中存在海洋洋流或是陸地水儲(chǔ)量的季節(jié)性周期信號(hào)殘留，對(duì)于第2、3、5、6模式，其對(duì)應(yīng)的時(shí)變模態(tài)主要是長(zhǎng)周期信號(hào)，信號(hào)的能量強(qiáng)度分別為13.23%、9.27%、3.29%、2.72%，其中可能包含有核幔邊界的特征信號(hào)。將重力殘項(xiàng)的時(shí)變模態(tài)與地磁場(chǎng)信號(hào)6個(gè)模式的時(shí)變模態(tài)進(jìn)行分析對(duì)比（圖14和圖17），結(jié)合有關(guān)學(xué)者的研究成果，認(rèn)為第3模式主要表征地核信號(hào)的變化特征［15-16］。

核幔邊界和地球表面地磁場(chǎng)第2模式的時(shí)變模態(tài)均與重力殘項(xiàng)第3模式的時(shí)變模態(tài)存在高度相關(guān)性，三者時(shí)變模態(tài)的信號(hào)特征如圖18所示。

圖18 重力殘項(xiàng)第3時(shí)變模態(tài)與地磁場(chǎng)（地球表面與核幔邊界）第2時(shí)變模態(tài)信號(hào)Fig.18 Third time-varying mode of gravity field and second time-varying mode of geomagnetic field (surface and core mantle boundary of the earth)

根據(jù)圖18可以看出，考慮其時(shí)間模態(tài)，對(duì)于重力殘項(xiàng)的第3模式與地磁場(chǎng)信號(hào)的第2模式而言，均存在較高的相關(guān)性；考慮其空間模態(tài)，結(jié)合圖5和圖8，地磁場(chǎng)徑向分量的二階導(dǎo)數(shù)在太平洋地區(qū)也存在較為明顯的變化信號(hào)，該變化特征在經(jīng)過(guò)PCA-ICA分解后得到凸顯（見(jiàn)圖13），尤其是地球表面的地磁場(chǎng)信號(hào)，其第1、2、4、6模式在太平洋地區(qū)均表現(xiàn)出較強(qiáng)的變化信號(hào)，與重力殘項(xiàng)的第3模式空間模態(tài)特征相對(duì)應(yīng)。

自2014年以來(lái)，太平洋地區(qū)地磁場(chǎng)徑向分量的二階導(dǎo)數(shù)存在著線性迅速變化的現(xiàn)象，局部地區(qū)的加速度變化在地球表面呈現(xiàn)東西向排列偶極子形式的簡(jiǎn)單結(jié)構(gòu)；對(duì)于重力信號(hào)的第3模式，在2014年以后，重力信號(hào)也呈現(xiàn)出線性迅速變化的現(xiàn)象，與磁場(chǎng)信號(hào)的時(shí)變特征存在著較高的相關(guān)性，但其第3空間模態(tài)在太平洋低緯度地區(qū)并未表現(xiàn)出明顯的“偶極子”形式的變化特征，其地區(qū)的空間模態(tài)存在較多的短波變化。

同時(shí)對(duì)比分析CHAOS-7模型重構(gòu)地球表面和核幔邊界的地磁場(chǎng)結(jié)果可知，目前將徑向分量二階導(dǎo)數(shù)從地球表面深入到核幔邊界仍存在一定難度，因?yàn)槠涔β首V在核幔邊界處隨著球諧階數(shù)的增加而增加，意味著占主導(dǎo)地位的短波長(zhǎng)成分并未受到良好的約束，此外，下地幔的電導(dǎo)率也存在不確定性，增加了向核幔邊界反演的難度，盡管如此，在太平洋中部低緯度地區(qū)地球表面的地磁場(chǎng)二階導(dǎo)數(shù)信號(hào)依然能在一定程度上反映核幔邊界相應(yīng)地區(qū)的一致性變化。

結(jié)合地磁場(chǎng)徑向分量二階導(dǎo)數(shù)信號(hào)（見(jiàn)圖5、圖8）及其第2模式的時(shí)變模態(tài)與空間模態(tài)（見(jiàn)圖18）可知，本實(shí)驗(yàn)成功探測(cè)出自2014年以后太平洋地區(qū)地磁場(chǎng)二階導(dǎo)數(shù)信號(hào)快速線性變化的特征，同時(shí)發(fā)現(xiàn)，重力殘項(xiàng)第3模式的時(shí)變模態(tài)的信號(hào)變化特征與地磁場(chǎng)徑向分量二階導(dǎo)數(shù)的時(shí)變模態(tài)的信號(hào)變化特征存在較好的相關(guān)性（見(jiàn)圖18）。

4 結(jié)語(yǔ)

基于2002年4月至2015年12月的GRACE衛(wèi)星數(shù)據(jù)，使用公開(kāi)發(fā)布的陸地水文模型（GLDAS模型）、海洋模型（ECCO模型）以及冰川均衡調(diào)整模型（ICE6G-D模型）對(duì)重力信號(hào)進(jìn)行改正，以消除陸地水文、陸地冰川、海洋質(zhì)量變化等非地核信號(hào)，從而分離出有效的地核信號(hào)。對(duì)于地磁場(chǎng)模型，采用CHAOS-7模型，反演出地球表面和核幔邊界的地磁地核場(chǎng)時(shí)變模型，時(shí)間跨度覆蓋2002年1月至2015年12月。采用PCA-ICA方法，對(duì)重力殘項(xiàng)信號(hào)與地磁場(chǎng)信號(hào)進(jìn)行統(tǒng)計(jì)分解，最后基于分解出的空間模態(tài)與時(shí)間模態(tài)對(duì)比分析兩者的相關(guān)性及一致性。

研究發(fā)現(xiàn)，重力殘項(xiàng)的第1、4模式為尚未扣除干凈的陸地水文信號(hào)和海洋洋流信號(hào)，存在明顯的周期性特征，其余模式存在長(zhǎng)周期特征，其中重力信號(hào)第3模式的時(shí)變特征與磁場(chǎng)信號(hào)第2模式的時(shí)變特征在同期內(nèi)存在相同的變化趨勢(shì)，兩者表現(xiàn)出較強(qiáng)的相關(guān)性，尤其在2014年以后，重力信號(hào)與磁場(chǎng)信號(hào)均呈現(xiàn)線性快速變化，與文獻(xiàn)［4-6］研究結(jié)果符合。磁場(chǎng)信號(hào)的空間模態(tài)反映了太平洋低緯度地區(qū)較為明顯的“偶極子”形式的變化特征。由于地磁地核場(chǎng)信號(hào)起源于地核，且實(shí)驗(yàn)所得重力殘項(xiàng)第3模態(tài)與地磁場(chǎng)信號(hào)表現(xiàn)出較強(qiáng)相關(guān)性，所以認(rèn)為本實(shí)驗(yàn)得到的重力信號(hào)中包含來(lái)自地核動(dòng)力學(xué)過(guò)程的特征信息。綜上所述，衛(wèi)星重力信號(hào)在經(jīng)過(guò)地表信號(hào)的剔除后，其包含的地核信號(hào)強(qiáng)度有所凸顯，也能在一定程度上反映核幔邊界處的動(dòng)力學(xué)過(guò)程。

作者貢獻(xiàn)聲明：

王正濤：提出研究主題、構(gòu)建研究框架、提出研究思路。

張勇剛：完善實(shí)驗(yàn)方法、收集實(shí)驗(yàn)數(shù)據(jù)、進(jìn)行數(shù)據(jù)實(shí)驗(yàn)、模型推演與論文撰寫(xiě)。