聯(lián)合時變重力數(shù)據(jù)與測高數(shù)據(jù)反演全球海平面變化及其分量貢獻(xiàn)

李 楊，郭金運，孫 玉，袁佳佳，常曉濤，張紅日

1. 山東科技大學(xué)測繪與空間信息學(xué)院，山東 青島 266590； 2. 福州大學(xué)空間數(shù)據(jù)挖掘與信息共享教育部重點實驗室，福建 福州 350108； 3. 自然資源部國土衛(wèi)星遙感應(yīng)用中心，北京 100048

全球變暖使地球氣候系統(tǒng)產(chǎn)生了一系列變化，例如地表氣溫升高、海洋熱含量增加、海洋浮冰融化和兩極冰蓋融化。在過去的50年間，溫室氣體排放使地球氣候系統(tǒng)中堆積了大量的熱量，海洋吸收了其中93%的熱量，剩余7%的熱量被土地、冰川、冰蓋和大氣吸收[1]。全球平均海平面(global mean sea level，GMSL)變化主要受海水體積熱膨脹、南極冰蓋融化、格陵蘭冰蓋融化、陸地冰川融化和陸地水儲量(terrestrial water storage，TWS)變化的影響。驗潮站數(shù)據(jù)顯示，19世紀(jì)中葉以后，全球平均海平面以1.2～1.9 mm/a的速率上升[2-4]。在GMSL上升和局部因素的作用下，沿海地區(qū)正受到海平面上升的威脅[5-7]。

目前，GMSL上升已成為地球氣候變化的一個重要響應(yīng)，研究GMSL有助于更好地觀測和解釋地球氣候變化，發(fā)現(xiàn)當(dāng)今地球觀測系統(tǒng)的不足[1]。政府間氣候變化專門委員會在第5次評估報告中指出[8]：1993年至2010年間使用衛(wèi)星測高數(shù)據(jù)估計的GMSL變化為3.2±0.4 mm/a，而使用其他數(shù)據(jù)估算的GMSL變化僅為2.8±0.5 mm/a，其中ANT、GRE和GLA的貢獻(xiàn)占比為50%，TWS減少導(dǎo)致的海平面上升占比為13%，海水溫鹽度變化引起的GMSL變化占比為37%，稱為比容(Steric)海平面變化。衛(wèi)星測高數(shù)據(jù)計算得到的海平面上升速率與其他數(shù)據(jù)估算的速率存在差異，這種差異可能反映出現(xiàn)今地球觀測系統(tǒng)不能完全捕捉到各因素對海平面變化的貢獻(xiàn)。

由于負(fù)荷自吸引效應(yīng)，地表或地球內(nèi)部的質(zhì)量遷移過程會引起大地水準(zhǔn)面變化和地表徑向位移，大地水準(zhǔn)面變化和地表位移呈現(xiàn)出的不均勻空間特征通常稱為海平面指紋(sea level fingerprint，SLF)或靜力平衡指紋[9-11]。隨著重力恢復(fù)與氣候試驗衛(wèi)星(gravity recovery and climate experiment，GRACE)的發(fā)射，SLF在海平面變化、地心運動和地球扁率變化等研究中得到了廣泛的應(yīng)用[12-14]。文獻(xiàn)[12,15]提出了一種聯(lián)合反演方法用于解釋全球海平面變化，該方法首先定義一組SLF和比容經(jīng)驗正交函數(shù)(empirical orthogonal function，EOF)，SLF和比容EOF分別表示不同的質(zhì)量遷移過程和海水比容對全球海平面變化造成的影響?；谶@組預(yù)定義的SLF和比容EOF，聯(lián)合GRACE數(shù)據(jù)和衛(wèi)星測高海平面異常(sea level anomaly，SLA)數(shù)據(jù)反演GMSL變化，并根據(jù)SLF和比容EOF將GMSL變化分解成ANT、GRE、GLA、TWS和Steric的貢獻(xiàn)。與傳統(tǒng)方法相比[16]，該反演方法可使用SLF將海水質(zhì)量變化引起的海平面變化進(jìn)一步細(xì)分成冰蓋融化、冰川融化、海陸水循環(huán)的貢獻(xiàn)；不過分依賴先驗的地球物理模型，可以同時估計冰川均衡調(diào)整(glacial isostatic adjustment，GIA)對絕對海平面變化的貢獻(xiàn)；無須預(yù)先扣除GRACE數(shù)據(jù)和衛(wèi)星測高數(shù)據(jù)中的GIA信號；不需要對GRACE數(shù)據(jù)進(jìn)行濾波，可以最大限度地保留GRACE數(shù)據(jù)的原始信號[12]。

文獻(xiàn)[12]實際反演得到的全球平均比容海平面(global mean steric sea level，GMSSL)變化僅為0.17±0.08 mm/a，其聯(lián)合反演結(jié)果中存在1.03±0.10 mm/a的殘差海平面變化無法解釋，反演效果不佳。文獻(xiàn)[15]利用有限元/體積海冰-海洋模型(finite-element/volume sea ice-ocean model，F(xiàn)ESOM)重新計算了一組比容EOF，對GMSL變化進(jìn)行了再次反演。其反演得的GMSSL變化為1.38±0.16 mm/a，該結(jié)果與直接使用海洋溫鹽數(shù)據(jù)產(chǎn)品計算的GMSSL變化相差較大，其估計的GMSSL變化占據(jù)GMSL變化的一半以上，目前普遍認(rèn)為海水質(zhì)量對GMSL變化的貢獻(xiàn)最大[15,17]。

本研究在文獻(xiàn)[12,15]的基礎(chǔ)上，使用EN4.2.1海洋溫鹽產(chǎn)品[18]計算了一組比容EOF，并根據(jù)相關(guān)輔助數(shù)據(jù)解算了一組SLF?；谶@組比容EOF和SLF，聯(lián)合GRACE數(shù)據(jù)、GRACE-FO數(shù)據(jù)和衛(wèi)星測高SLA數(shù)據(jù)，反演了2002年4月至2020年2月的全球海平面變化，將海平面變化分解成ANT、GRE、GLA、TWS、GIA、Steric的貢獻(xiàn)。以期解決文獻(xiàn)[12,15]中估算的GMSSL變化與前人相關(guān)文獻(xiàn)不一致的問題，使用聯(lián)合反演方法獲得合理的結(jié)果。

1 數(shù)據(jù)與方法

1.1 SLF和比容EOF

本文定義SLF的方式與文獻(xiàn)[15]定義SLF的方式一致，假設(shè)某個地理區(qū)域發(fā)生質(zhì)量變化，其造成的大地水準(zhǔn)面變化、地表徑向位移和相對海平面變化稱為SLF。本研究使用的SLF用于區(qū)分不同的質(zhì)量遷移過程，使用所有的SLF同時擬合GRACE時變重力數(shù)據(jù)或GRACE后繼衛(wèi)星的(GRACE Follow-On，GRACE-FO)時變重力數(shù)據(jù)可以得到每個SLF在時間域的變化。SLF不在于數(shù)量的多少，而在于使用的SLF能否恢復(fù)重力場的變化。

本文將地表質(zhì)量遷移過程分為兩類，一類是冰蓋和冰川的融化，另一類是TWS變化。研究中把第一類質(zhì)量遷移過程所在的位置劃分成一些地理區(qū)域，假設(shè)每個地理區(qū)域發(fā)生單位質(zhì)量變化，然后對每個地理區(qū)域求解一個SLF。根據(jù)文獻(xiàn)[19—20]對南極冰蓋匯水區(qū)的定義，將南極冰蓋劃分成了25個地理區(qū)域；根據(jù)格陵蘭冰蓋匯水區(qū)的定義并以2000 m的高程為界限，將格陵蘭冰蓋劃分成15個地理區(qū)域；陸地冰川的位置由全球陸地冰雪空間觀測項目提供[21]，劃分了14個陸地冰川區(qū)域。TWS變化較為復(fù)雜，難以手動劃分地理區(qū)域，為了計算TWS相關(guān)的SLF，本文對Water GAP全球水文模型(water GAP global hydrology model，WGHM)進(jìn)行經(jīng)驗正交函數(shù)分析提取了前60個EOF[22-23]，這些EOF的累計方差貢獻(xiàn)率在90%以上。在提取EOF之前剔除了WGHM水文模型中的冰川區(qū)域，對每個EOF求解了一個SLF。水文模型不能捕捉到人為因素引起的TWS變化，為了彌補(bǔ)水文模型和GRACE/GRACE-FO數(shù)據(jù)的不一致，本文額外引入了3個SLF捕捉里海、北印度平原和維多利亞湖附近區(qū)域的地下水減少信號，一共解算了63個TWS相關(guān)的SLF。此外，研究中額外定義了一組SLF以考慮GIA對海平面變化的影響：基于ICE-6G冰川歷史模型[24]，解算了2個Fennoscandia區(qū)域的SLF，3個Laurentide區(qū)域的SLF；基于南極的IJ05冰川歷史模型[25]，解算了1個Antarctica的SLF。使用這6個SLF共同表示GIA對全球大地水準(zhǔn)面變化和地表徑向位移的影響。

本文使用英國氣象局(Met Office)發(fā)布的EN4.2.1格網(wǎng)化海洋溫鹽產(chǎn)品并借助開源程序Gibbs-SeaWater海洋學(xué)工具箱計算比容海平面變化[18]，然后對比容海平面變化進(jìn)行經(jīng)驗正交函數(shù)分析提取比容EOF以參數(shù)化海水比容變化對全球海平面的影響。EN4.2.1海洋溫鹽產(chǎn)品融合了全球海洋數(shù)據(jù)庫13(world ocean database 13，WOD13)、全球溫度和鹽度剖面項目(global temperature and salinity profile project，GTSPP)、Argo計劃和北極概要流域觀測項目(arctic synoptic basin-wide observations，ASBO)，提供了1900年至今的海洋三維溫鹽場信息，相對于僅包含Argo數(shù)據(jù)的溫鹽產(chǎn)品，EN4.2.1海洋溫鹽產(chǎn)品包含的比容海平面變化信息更多。研究中使用了1990年1月至2020年12月的EN4.2.1海洋溫鹽產(chǎn)品，EN4.2.1海洋溫鹽產(chǎn)品并未用作觀測數(shù)據(jù)，而是用于提取比容EOF，一共提取了239個比容EOF，其累計方差貢獻(xiàn)率在99%以上。使用1990年至2020年的比容海平面變化提取比容EOF能使比容EOF包含的比容海平面變化空間特征更多，例如1997年至1998年的強(qiáng)El Nio事件對比容海平面變化的影響也被囊括在內(nèi)。

在產(chǎn)業(yè)發(fā)展初期，太陽能電池組件成本和發(fā)電價格是光伏發(fā)電應(yīng)用和推廣的障礙，主要依賴政府持續(xù)穩(wěn)定的扶持政策才能發(fā)展壯大。各國紛紛出臺“可再生能源法”采用法制手段強(qiáng)制太陽能光伏并網(wǎng)發(fā)電，為光伏產(chǎn)業(yè)快速發(fā)展奠定法律基礎(chǔ)；通過政策扶持使生產(chǎn)者、供應(yīng)者和使用者都得到很多的補(bǔ)貼，有利于推進(jìn)本國光伏產(chǎn)業(yè)的發(fā)展；鼓勵推進(jìn)屋頂光伏發(fā)電工程，促進(jìn)整個光伏產(chǎn)業(yè)迅速發(fā)展。

1.2 衛(wèi)星觀測數(shù)據(jù)

采用美國德克薩斯大學(xué)空間研究中心(Center for Space Research，CSR)發(fā)布的GRACE/GRACE-FO RL06版本的月時變重力場模型，該模型為截斷至60階次的球諧系數(shù)產(chǎn)品，已經(jīng)扣除了非潮汐大氣和海洋信號[26]。研究中使用了2002年4月至2020年2月共182個月的時變重力場模型，未替換GRACE/GRACE-FO數(shù)據(jù)中缺失的C10、C11、S11系數(shù)和不準(zhǔn)確的C20系數(shù)，已有研究表明基于海平面指紋的反演方法可以估計這些系數(shù)的變化[13-14]。2016年8月以后，受衛(wèi)星單加速度計工作模式影響，GRACE/GRACE-FO數(shù)據(jù)的C30系數(shù)不準(zhǔn)確，得益于2012年激光相對論衛(wèi)星(laser relativity satellite，LARES)的發(fā)射，衛(wèi)星激光測距(satellite laser ranging，SLR)數(shù)據(jù)提供了高精度的C30系數(shù)，本文采用戈達(dá)德航天中心(Goddard Space Flight Center，GSFC)解算的SLRC30系數(shù)替換了2016年8月以后GRACE/GRACE-FO產(chǎn)品的C30系數(shù)[27]，該C30系數(shù)與原始GRACE/GRACE-FO的C30系數(shù)一致性好，已經(jīng)被CSR用于改正其生產(chǎn)的Mascon產(chǎn)品。根據(jù)這182個月的時變重力場模型計算了研究時間段的平均重力場，將每個月的時變重力場模型都扣除這個平均重力場再乘以平均地球半徑以表示大地水準(zhǔn)面變化，作為聯(lián)合反演的輸入數(shù)據(jù)。

本文采用AVISO(archiving,validation and interpretation of satellite oceanographic data)提供的2002年4月至2020年2月的格網(wǎng)化月均SLA數(shù)據(jù)。該數(shù)據(jù)融合了多顆衛(wèi)星測高數(shù)據(jù)，進(jìn)行了必要的地球物理改正及大氣改正，包括電離層改正、干濕對流層校正、固體潮、海洋潮汐、海況偏差和逆氣壓改正等，該數(shù)據(jù)的格網(wǎng)大小為0.25°×0.25°。由于時變重力數(shù)據(jù)在某些月份缺失，為了保證聯(lián)合反演中時變重力數(shù)據(jù)和衛(wèi)星測高數(shù)據(jù)的參考時間一致，本研究剔除了缺失時變重力數(shù)據(jù)當(dāng)月的格網(wǎng)SLA數(shù)據(jù)，然后計算研究時間段內(nèi)SLA數(shù)據(jù)的均值，并從每月的SLA數(shù)據(jù)中扣除該均值。

本文沒有恢復(fù)GRACE/GRACE-FO時變重力數(shù)據(jù)的非潮汐大氣和海洋信號，如果恢復(fù)該信號，那么在聯(lián)合反演中，需要額外引入一組基函數(shù)表示非潮汐大氣和海洋信號。為保持變重力數(shù)據(jù)和衛(wèi)星測高數(shù)據(jù)的信號一致，使用AOD1B(atmosphere and ocean de-aliasing level-1B)模型扣除了SLA數(shù)據(jù)中的非潮汐海洋信號[26,28]?？鄢鼳OD1B模型之前，每個月的AOD1B模型需扣除其在海洋上的空間均值。非潮汐海洋信號對GMSL變化沒有貢獻(xiàn)，在研究GMSL變化時可以不恢復(fù)該信號，但是該信號對區(qū)域海平面變化有影響，在研究區(qū)域海平面變化時，必須恢復(fù)該信號。考慮到時變重力數(shù)據(jù)和衛(wèi)星測高數(shù)據(jù)的空間分辨率差異，本研究將SLA數(shù)據(jù)重采樣成1°×1°的格網(wǎng)數(shù)據(jù)，作為聯(lián)合反演的輸入數(shù)據(jù)。

1.3 聯(lián)合反演方法

本文所使用的方法稱為聯(lián)合反演方法，其基本假設(shè)是造成全球海平面變化的信號在空間上可以較好地分離并建模，并且各類信號造成的空間特征彼此線性無關(guān)，不同信號造成的空間特征由SLF或者比容EOF表示。

時變重力數(shù)據(jù)表示的大地水準(zhǔn)面變化是由質(zhì)量遷移過程導(dǎo)致的，本文認(rèn)為大地水準(zhǔn)面變化是SLF的線性組合，每月可以在頻域建立一組時變重力觀測方程

Gx(t)SLF=b(t)grace+ε

(1)

式中，b(t)grace是GRACE/GRACE-FO數(shù)據(jù)表示的大地水準(zhǔn)面變化；設(shè)計矩陣G由SLF組成，每個列向量是一個SLF，表示相應(yīng)的質(zhì)量遷移過程對大地水準(zhǔn)面造成的影響，設(shè)計矩陣G和觀測向量b(t)grace中的元素均為球諧系數(shù)且相同行系數(shù)的階次一致；解向量x(t)SLF表示SLF在時間域的變化。

本文沒有填補(bǔ)GRACE/GRACE-FO數(shù)據(jù)的C10、C11和S11系數(shù)，而是剔除了觀測方程中的C10、C11和S11系數(shù)，在聯(lián)合反演中估計它們的變化[14]?？紤]到GRACE/GRACE-FO數(shù)據(jù)的C20系數(shù)不準(zhǔn)確，還額外剔除了觀測方程中的C20系數(shù)，C20系數(shù)也可以基于海平面指紋反演得到[13]。嘗試了利用聯(lián)合反演方法估計C30系數(shù)的變化，但結(jié)果并不理想，因此在數(shù)據(jù)預(yù)處理階段使用了SLR的C30系數(shù)代替了2016年8月以后GRACE/GRACE-FO的C30系數(shù)。

絕對海平面會受質(zhì)量遷移引起的大地水準(zhǔn)面變化而變化，還會受到海水溫鹽度、海洋環(huán)流、大氣等影響。在文中非潮汐海洋信號已經(jīng)從SLA數(shù)據(jù)中扣除，因此，本文認(rèn)為絕對海平面變化可由預(yù)定義的SLF和比容EOF線性表示，每月可在空間域建立一組衛(wèi)星測高觀測方程

Ax(t)=b(t)altimetry+ε

(2)

式中，觀測向量b(t)altimetry是SLA數(shù)據(jù)表示的絕對海平面變化；矩陣B由SLF組成，每個列向量是一個SLF，矩陣C由比容EOF組成，每個列向量是一個比容EOF，設(shè)計矩陣A由矩陣B和矩陣C組成；x(t)EOF表示比容EOF在時間域的變化，解向量x(t)由x(t)SLF和x(t)EOF組成，表示所有基函數(shù)在時間域的變化。該觀測方程中使用的SLF和比容EOF均由格網(wǎng)數(shù)據(jù)表示，矩陣B、矩陣C和觀測向量b(t)altimetry中相同行元素的經(jīng)緯度一致。對于格網(wǎng)大小為1°×1°的SLA數(shù)據(jù)，觀測向量b(t)altimetry有64 800個元素，需要注意的是，衛(wèi)星測高觀測方程中應(yīng)僅含有海洋上的有效格網(wǎng)點值，本研究剔除了陸地上的格網(wǎng)點和格網(wǎng)值為空值的格網(wǎng)點，未剔除近海區(qū)域的測高數(shù)據(jù)。剔除無效數(shù)據(jù)之后，每月的觀測向量b(t)altimetry約有33 000個元素。

(3)

式中

(4)

ATPaltimetryb(t)altimetry

(5)

2 結(jié)果與討論

2.1 反演結(jié)果分析

根據(jù)聯(lián)合反演的結(jié)果重建了全球平均絕對海平面(global mean absolute sea level，GMASL)變化。為驗證聯(lián)合反演的結(jié)果，本文還根據(jù)觀測向量b(t)altimetry計算了GMASL變化，并根據(jù)聯(lián)合反演中觀測向量b(t)altimetry的殘差計算了全球平均的殘差時間序列，如圖1所示。GRACE/GRACE-FO數(shù)據(jù)在某些月份缺失，無法進(jìn)行聯(lián)合反演，因此本文展示的時間序列中含有缺失值。這些時間序列在2017年7月—2018年5月連續(xù)斷開11個月，本文未對該時段的缺失值做處理，僅對其他的缺失值進(jìn)行了線性插值補(bǔ)全。

圖1 全球平均絕對海平面變化Fig.1 Global mean absolute sea level change

由圖1可知，殘差時間序列具有明顯的周期特征，擬合得知其周年振幅為1.29±0.21 mm。聯(lián)合反演和衛(wèi)星測高數(shù)據(jù)得到的GMASL變化的線性趨勢分別為3.01±0.09 mm/a和3.29±0.10 mm/a，殘差序列的線性趨勢為0.28±0.02 mm/a。從線性趨勢看，當(dāng)前SLF和比容EOF配置下的聯(lián)合反演可以恢復(fù)超過90%的GMASL變化，僅有一小部分海平面變化無法解釋。

聯(lián)合反演中觀測向量b(t)altimetry殘差的標(biāo)準(zhǔn)差如圖2(a)所示。由圖可知，聯(lián)合反演可以重建大空間尺度的海平面變化，但是在一些海域反演的效果不理想。研究發(fā)現(xiàn)聯(lián)合反演殘差的標(biāo)準(zhǔn)差與SLA數(shù)據(jù)的標(biāo)準(zhǔn)差在空間上存在相關(guān)性，如圖2(b)所示。這表明聯(lián)合反演的效果受SLA數(shù)據(jù)的可變性影響，在SLA變化大的海域，聯(lián)合反演的效果不好，SLA的高可變性主要與洋流和洋流導(dǎo)致的渦流等海洋現(xiàn)象有關(guān)。受日本暖流、千島寒流和北太平洋暖流的影響，聯(lián)合反演在日本東部海域的反演效果較差；受墨西哥灣暖流的影響，聯(lián)合反演在美國東部海域的反演效果較差；受東澳大利亞暖流的影響，聯(lián)合反演在澳大利亞東部海域的反演效果較差。衛(wèi)星測高數(shù)據(jù)在數(shù)據(jù)預(yù)處理階段使用AOD1B模型扣除了非潮汐海洋信號，由于AOD1B模型本身存在誤差且其空間分辨率較低，這樣的操作可能影響反演的結(jié)果。如果在聯(lián)合反演中加回AOD1B模型，對AOD1B模型或其他海洋模型提取EOF以表示非潮汐海洋信號，同時估計其時變，有望能減少海洋現(xiàn)象顯著海域的反演殘差。

圖2 聯(lián)合反演中殘差標(biāo)準(zhǔn)差的空間分布與衛(wèi)星測高數(shù)據(jù)標(biāo)準(zhǔn)差的空間分布Fig.2 The spatial distribution of standard deviation of residuals in joint inversion and the spatial distribution of standard deviation of satellite altimetry data

為進(jìn)一步驗證方法的可靠性，本文將聯(lián)合反演得到的C10、C11、S11和C20系數(shù)與僅使用式(1)進(jìn)行反演得到的結(jié)果[13]、GSFC發(fā)布的SLR結(jié)果[27]和文獻(xiàn)[30]得到的結(jié)果(GRACE-OBP方法)進(jìn)行了對比。C10、C11和S11系數(shù)的時間序列如圖3所示，從上到下依次為C10、C11和S11，各時間序列均未加回AOD1B模型。C20系數(shù)的時間序列如圖4所示，虛線表示相應(yīng)的時間序列進(jìn)行12個月的滑動平均濾波得到的結(jié)果，對時間序列進(jìn)行了上下平移以便辨識，各時間序列均未加回AOD1B模型。

由圖3可知，聯(lián)合反演得到的3個一階項系數(shù)與GRACE-OBP方法的一階項系數(shù)在周年振幅方面比較一致，3種方法的C10和C11系數(shù)的線性趨勢基本一致，基于海平面指紋反演得到的S11系數(shù)的線性趨勢均大于GRACE-OBP的結(jié)果。C10、C11和S11系數(shù)的線性趨勢、周年振幅和周年相位見表1。由圖4可知，聯(lián)合反演、GRACE-OBP方法和僅使用GRACE/GRACE-FO進(jìn)行指紋法反演得到的C20系數(shù)較為一致，SLR的C20時間序列振幅小于其他3種時間序列的振幅，年際變化基本一致。C20時間序列的線性趨勢、周年振幅和周年相位，見表2。由擬合結(jié)果可知，聯(lián)合反演得到的C20時間序列在線性趨勢和周年相位兩個方面與GRACE-OBP方法和SLR的結(jié)果一致，SLR的C20周年振幅小于聯(lián)合反演和GRACE-OBP方法得到的周年振幅。

圖3 使用不同方法得到的C10、C11和S11時間序列Fig.3 Time series of C10, C11 and S11 derived from different approaches

圖4 C20時間序列Fig.4 Time series of C20

表1 C10、C11和S11系數(shù)的線性趨勢、周年振幅和周年相位Tab.1 Linear trend, annual amplitude and annual phase of time series of C10, C11 and S11

表2 C20時間序列的線性趨勢、周年振幅和周年相位Tab.2 Linear trend, annual amplitude and annual phase of time series of C20

2.2 全球平均海平面變化

根據(jù)聯(lián)合反演的結(jié)果計算了ANT、GRE、GLA、TWS、Steric對GMSL和GMASL的貢獻(xiàn)。GMSL變化及其分量的貢獻(xiàn)如圖5所示，劃分成6個分量的貢獻(xiàn)，擬合并去除了時間序列中的年周期項和半年周期項，進(jìn)行了6個月的滑動平均濾波，對時間序列進(jìn)行了上下平移以便辨識。各時間序列的線性趨勢、周年振幅和周年相位，見表3。由圖5可知，冰蓋、冰川融化對GMSL變化的貢獻(xiàn)較為穩(wěn)定，TWS和Steric對GMSL變化的貢獻(xiàn)波動較大。2004—2014年TWS時間序列的線性趨勢為0.05±0.08 mm/a，這表明TWS對GMSL變化幾乎沒有長期貢獻(xiàn)，2014年后TWS時間序列出現(xiàn)上升的轉(zhuǎn)折點，這暗示2014年后全球陸地水儲量的減少。據(jù)表3可知，冰蓋和冰川融化是當(dāng)今GMSL上升的主要原因，它們對GMSL貢獻(xiàn)的線性趨勢達(dá)到1.68 mm/a，受全球氣候變暖影響，海水溫鹽度變化導(dǎo)致的海水體積膨脹是GMSL上升的次要原因。比容作用和TWS變化是主導(dǎo)GMSL周期性變化的主要因素，其周年振幅分別達(dá)到了8.66±0.34 mm和4.01±0.44 mm。殘差時間序列的相位為61±3 d，與GMSSL的相位(75±4 d)接近，海水質(zhì)量變化在聯(lián)合反演中被GRACE/GRACE-FO數(shù)據(jù)約束，從這兩方面看殘差屬于比容信號的可能性較大。

表3 擬合得到的線性趨勢、周年振幅和周年相位Tab.3 Estimated annual amplitudes, phases and linear trends

圖5 全球平均海平面變化Fig.5 Global mean sea level change

本文除了對1990—2020年的比容海平面變化進(jìn)行經(jīng)驗正交函數(shù)分析提取比容EOF外，還對2000—2020年的比容海平面變化提取了比容EOF，分析了不同比容EOF對聯(lián)合反演的影響，兩組比容EOF的累計方差貢獻(xiàn)率均為99%。2條反演得到的GMSSL時間序列與使用EN4.2.1海洋溫鹽產(chǎn)品計算的GMSSL時間序列如圖6所示。

圖6 聯(lián)合反演和EN4.2.1海洋溫鹽產(chǎn)品計算得到的全球平均比容海平面變化Fig.6 Global mean steric sea level change derived from joint inversion and EN4.2.1 product

由圖6可知，相較于使用EN4.2.1海洋溫鹽產(chǎn)品計算得的GMSSL變化，利用不同比容EOF進(jìn)行聯(lián)合反演得到的GMSSL在年際變化和季節(jié)性變化方面更加接近。EN4.2.1海洋溫鹽產(chǎn)品計算得到的GMSSL在2002—2006年的年際變化明顯，而同期聯(lián)合反演解算得到的GMSSL不存在明顯的年際變化，Argo計劃早期浮標(biāo)剖面數(shù)據(jù)少，聯(lián)合反演方法對GMSSL的估計可能更加合理。

擬合得到了圖6中未經(jīng)濾波時間序列的線性趨勢、周年振幅和周年相位，結(jié)果見表4，兩組聯(lián)合反演得到的GMSSL變化的周年振幅均大于使用EN4.2.1海洋溫鹽產(chǎn)品計算得到的振幅。受淺海溫鹽剖面數(shù)據(jù)少、缺少深層海水溫鹽信息等缺點，實測海洋溫鹽數(shù)據(jù)可能低估比容海平面變化[31]，本文認(rèn)為聯(lián)合反演的比容EOF具備捕捉深海的溫鹽度變化信息的能力，也能彌補(bǔ)早期剖面數(shù)據(jù)不足的缺點，因此得到的GMSSL變化應(yīng)大于直接計算的結(jié)果。然而，使用2000—2020年的比容EOF進(jìn)行聯(lián)合反演得到的GMSSL變化的線性趨勢僅為0.87±0.04 mm/a，小于使用EN4.2.1海洋溫鹽產(chǎn)品計算的結(jié)果，比容EOF并沒有從衛(wèi)星測高數(shù)據(jù)中捕捉到足夠的比容海平面變化信息，這說明2000—2020年的比容EOF包含的比容信號過少，使用1990—2020年的比容EOF得到的結(jié)果更合理一些。

表4 全球平均比容海平面變化的線性趨勢、周年振幅和周年相位Tab.4 Linear trend, annual amplitude and annual phase of global mean steric sea level change

使用1990—2020年的比容EOF反演得到的GMSSL變化的線性趨勢為1.08±0.05 mm/a，大約占GMSL變化的三分之一，相比于文獻(xiàn)[12]和文獻(xiàn)[15]研究估計的GMSSL，本文結(jié)果與相關(guān)文獻(xiàn)的記載的結(jié)果更符合[32-33]。聯(lián)合反演方法對GMSSL變化的估計很大程度上取決于先驗比容EOF包含比容信號的多少。文獻(xiàn)[12]首次反演得到的GMSSL變化僅為0.17±0.08 mm/a，本文認(rèn)為是因為其使用ishii海洋溫鹽產(chǎn)品計算0～700 m的比容海平面變化導(dǎo)致的，這使得比容EOF包含的比容變化信號過少，直接導(dǎo)致了比容EOF無法從測高數(shù)據(jù)中捕捉到完整的比容海平面變化。文獻(xiàn)[15]使用FESOM模型提取了新的比容EOF用于反演，得到的GMSSL變化為1.38±0.16 mm/a，超過了海水質(zhì)量變化的貢獻(xiàn)(1.06 mm/a)，本文認(rèn)為這是由于其解算的比容EOF包含一部分噪聲信息，噪聲信息可能與SLF存在一定的相關(guān)性，導(dǎo)致反演時錯誤地估計了海水比容效應(yīng)和海水質(zhì)量變化。

在GMSL變化研究中通常使用先驗的GIA模型扣除GRACE數(shù)據(jù)和衛(wèi)星測高數(shù)據(jù)中的GIA信號，或扣除一個大約-0.3 mm/a的GIA改正值將GMASL轉(zhuǎn)換成GMSL[34-35]。本聯(lián)合反演得到的GIA對GMASL的貢獻(xiàn)僅為-0.06 mm/a，一方面這是因為研究中使用的南極地區(qū)的GIA SLF是基于IJ-05冰川歷史模型解算得到，相比于ICE系列模型，IJ-05模型估計的南極GIA信號對現(xiàn)今GMASL貢獻(xiàn)的量級較低[36]，另一方面，是因為本研究解算得到的南極區(qū)域GIA信號量級比先驗的IJ-05模型量級小。這一結(jié)果與文獻(xiàn)[15]的研究類似，其估計的GIA貢獻(xiàn)為-0.1 mm/a，雖然使用的SLF不同，但是在聯(lián)合反演中解算得到的南極GIA信號都較低。

利用聯(lián)合反演結(jié)果計算了全球平均海水質(zhì)量(global mean ocean mass，GMOM)變化，作為對照，本文依據(jù)式(1)僅使用GRACE/GRACE-FO數(shù)據(jù)反演了GMOM變化，還間接地基于測高數(shù)據(jù)和海洋溫鹽數(shù)據(jù)計算了GMOM變化，結(jié)果如圖7所示。2006—2016年間，測高數(shù)據(jù)減比容數(shù)據(jù)估算的GMOM變化與反演得到的GMOM變化基本一致，殘差都較小。2016年8月后，受衛(wèi)星單加速度計工作模式的影響，僅使用GRACE/GRACE-FO時變重力數(shù)據(jù)的反演可能低估了GMOM變化，這一現(xiàn)象與文獻(xiàn)[37]吻合。值得注意的是，融合了衛(wèi)星測高數(shù)據(jù)的聯(lián)合反演在一定程度上減弱了GMOM變化被低估的現(xiàn)象。在GRACE-FO數(shù)據(jù)時期，將僅使用時變重力數(shù)據(jù)反演的GMOM變化與基于測高數(shù)據(jù)和比容數(shù)據(jù)估算的GMOM變化相減，其殘差的RMS為13.60 mm，而聯(lián)合反演估計的GMOM變化與基于測高數(shù)據(jù)和比容數(shù)據(jù)估算的GMOM變化相減，其殘差的RMS僅為5.78 mm。

圖7 全球平均海水質(zhì)量變化(以等效海平面變化表示)Fig.7 Global mean ocean mass change(expressed as equivalent sea level change)

2.3 區(qū)域海平面變化

本文使用的聯(lián)合反演方法理論上可以分解海洋上任意位置的海平面變化，上文已經(jīng)展示了在近海區(qū)域聯(lián)合反演殘差的標(biāo)準(zhǔn)差沒有明顯增加的特點，這表明該方法不僅適用于GMSL變化的研究，還可用于區(qū)域海平面變化的研究[38]。本文試驗恢復(fù)了非潮汐海洋信號，選取了20個海域計算區(qū)域的相對海平面變化，并將海平面變化分解成ANT、GRE、GLA、TWS、GIA、Steric、Other的貢獻(xiàn)，Other指的是聯(lián)合反演的殘差和AOD1B模型的貢獻(xiàn)，結(jié)果如圖8所示。

圖8 各區(qū)域的相對海平面變化Fig.8 Relative sea level change in some regions

不同于GMSL變化，區(qū)域相對海平面變化還會受GIA影響，近海的相對海平面變化通常還受地質(zhì)構(gòu)造、地震、陸地沉降等影響[39]。本文并未顧及這些區(qū)域信號，但可在聯(lián)合反演的后處理中加上這些信號。由圖8可知，雖然海水比容效應(yīng)是GMSL變化的次要因素，但在很多近海區(qū)域，如墨西哥西部海域、印度西部海域、塔斯曼海等，比容變化的貢獻(xiàn)更大，超過了海水質(zhì)量變化的貢獻(xiàn)。

本文得到的中國南海相對海平面變化為4.33±1.81 mm/a，其中海水比容變化的貢獻(xiàn)占1.38 mm/a，ANT、GRE和GLA冰融化的貢獻(xiàn)分別為0.47、0.79和0.69 mm/a，TWS的貢獻(xiàn)為0.34 mm/a。中國東海區(qū)域的反演效果不佳，這可能與中國近海Argo剖面數(shù)據(jù)較少有關(guān)。在北海、歐洲西部海域和地中海聯(lián)合反演的質(zhì)量最差，殘差變化的線性趨勢分別為2.48、2.12和1.40 mm/a，澳大利亞東部海域和塔斯曼海的反演效果最好，其殘差趨勢僅為0.02 mm/a和0.12 mm/a。由圖8中的20個海域反演情況來看，僅有歐洲的3個海域和中國東海反演效果較差，相比于GMSL變化0.28 mm/a的殘差，聯(lián)合反演中其他16個海域的Other貢獻(xiàn)在可接受的范圍內(nèi)，表明本研究結(jié)果可以得到對GMSL變化和區(qū)域海平面變化的一致性估計，可適用于近海的海平面變化研究。

3 結(jié) 論

本文基于一組預(yù)定義的SLF和比容EOF，聯(lián)合2002年4月—2020年2月的時變重力數(shù)據(jù)和衛(wèi)星測高數(shù)據(jù)反演了全球海平面變化，在最小二乘意義上給出了ANT、GRE、GLA、TWS、GIA和Steric對GMSL和GMASL變化的貢獻(xiàn)。

結(jié)果顯示，本文解算的SLF和比容EOF可以解釋90%以上的GMASL，2002—2020年GMASL和GMSL的線性趨勢分別為3.29±0.09 mm/a和3.47±0.09 mm/a。聯(lián)合反演得到的GMSSL的線性趨勢為1.08±0.05 mm/a，與相關(guān)文獻(xiàn)中的結(jié)果一致，比文獻(xiàn)[12]和文獻(xiàn)[15]估計的GMSSL變化合理。研究發(fā)現(xiàn)，結(jié)合衛(wèi)星測高數(shù)據(jù)和時變重力數(shù)據(jù)的聯(lián)合反演方法能夠減弱GRACE-FO衛(wèi)星時期海水質(zhì)量變化被低估的現(xiàn)象。本文分析了聯(lián)合反演中衛(wèi)星測高觀測方程中觀測向量殘差的標(biāo)準(zhǔn)差，發(fā)現(xiàn)聯(lián)合反演在海洋動力現(xiàn)象顯著的海域效果較差，在近海區(qū)域反演效果并沒有降低?；诖税l(fā)現(xiàn)，利用該聯(lián)合反演方法分解了20個區(qū)域的相對海平面變化，將該聯(lián)合反演方法用于區(qū)域海平面變化的研究。