基于Slepian方法和地面重力觀測確定時變重力場模型：以2011—2013年華北地區(qū)數(shù)據(jù)為例

韓建成， 陳石* ， 盧紅艷， 徐偉民

1 中國地震局地球物理研究所, 北京 100081 2 北京白家疃地球科學(xué)國家野外科學(xué)觀測研究站, 北京 100095

0 引言

時變重力場反映地球系統(tǒng)內(nèi)物質(zhì)的分布、運(yùn)動和變化，是研究地球系統(tǒng)內(nèi)物質(zhì)運(yùn)動狀態(tài)及其動力學(xué)機(jī)制的重要約束信息，可為解決人類面臨的資源、環(huán)境和災(zāi)害等緊迫問題提供重要的基礎(chǔ)地球物理信息(許厚澤等, 2012; 寧津生等, 2016).近十幾年來，GRACE(Gravity Recovery And Climate Experiment)衛(wèi)星已成功獲取全球均勻覆蓋的時變重力觀測數(shù)據(jù)，使地球時變重力場模型在分辨率和精度上實(shí)現(xiàn)了重大提升(Tapley et al., 2004; Flechtner et al., 2014; 寧津生等, 2016).GRACE月時變重力場在大空間尺度上定量揭示了全球環(huán)境變化，例如：陸地水儲量變化(Rodell et al., 2018)、冰川消融(Harig and Simons, 2015; 高春春等, 2019)、海洋變化(Reager et al., 2016; Jin et al., 2020)、冰后回彈(Riva et al., 2009)以及強(qiáng)地震(Han et al., 2006; Panet et al., 2018)等導(dǎo)致的地球表層質(zhì)量重新分布.

GRACE重力衛(wèi)星具有可重復(fù)、大范圍的觀測優(yōu)勢，但其軌道距地面約500 km，所觀測到的地球重力場信號已大幅衰減，其中短波長(高頻)部分衰減幅度最大(周新等, 2011; 徐新禹等, 2018).因此，GRACE對重力場信號的中長波長部分較為敏感，對短波長部分敏感性較弱，不利于恢復(fù)較高階的時變地球重力場模型.GRACE的實(shí)際空間分辨率約400 km (Tapley et al., 2004)，其數(shù)據(jù)存在空間域的南北向條帶等噪聲需做濾波處理(Swenson and Wahr, 2006; Chen et al., 2007)，會進(jìn)一步降低空間分辨率，難以滿足環(huán)境保護(hù)、資源監(jiān)測以及災(zāi)害預(yù)測等研究對高精度高分辨率時變重力場的需求(許厚澤等, 2012; Alley and Konikow, 2015; Feng et al., 2018).

地面重力時變觀測在地表進(jìn)行，通過定期同點(diǎn)位復(fù)測來獲得重力場隨時間的變化.與GRACE重力衛(wèi)星相比，地面重力測量的優(yōu)點(diǎn)在于重力信號衰減幅度小，觀測結(jié)果的短波長部分更加精確(Bomfim et al., 2013)，空間分辨率也更高(Van Camp et al., 2017)，適于分析區(qū)域性、近地表的物質(zhì)運(yùn)動和遷移現(xiàn)象(許才軍和尹智, 2014; Van Camp et al., 2017).但地面重力觀測易受地理和環(huán)境等因素制約：一方面，很多區(qū)域(如南極、沙漠和冰川等)難以施測，存在數(shù)據(jù)空白區(qū)，無法像GRACE衛(wèi)星觀測那樣實(shí)現(xiàn)數(shù)據(jù)全球覆蓋；另一方面，在可施測區(qū)域內(nèi)，測點(diǎn)空間位置也難以實(shí)現(xiàn)均勻分布.直接利用地面重力觀測數(shù)據(jù)確定區(qū)域重力場模型時，由于數(shù)據(jù)僅區(qū)域覆蓋且空間不均勻分布，無法滿足球諧分析這樣的經(jīng)典重力場恢復(fù)理論對數(shù)據(jù)的要求(基函數(shù)在區(qū)域內(nèi)部不滿足正交性).因此，需借助局部重力場恢復(fù)理論和方法.

目前應(yīng)用較為廣泛的局部重力場恢復(fù)方法包括最小二乘配置法(Hwang et al., 2014; 王灼華等, 2017)、球冠諧分析法(Li et al., 1995; 王燚和姜效典, 2017)、Mascon (Mass concentration)法(Watkins et al., 2015; 鄒賢才等, 2016; Ran et al., 2018)等.除上述方法外，國內(nèi)外學(xué)者還發(fā)展了多種采用特殊基函數(shù)的區(qū)域重力場建模方法，如矩諧分析法、徑向基函數(shù)(radial basis function)法以及 Slepian 局部譜分析法等.矩諧分析方法選取三角函數(shù)作為基函數(shù)，便于計(jì)算和展開到較高階次(蔣濤等, 2014).徑向基函數(shù)方法可定義多種基函數(shù)，具備在空間域和頻率域局部化的特性(Klees et al., 2008; 吳懌昊等, 2016; 楊帆等, 2017).Slepian 局部譜分析方法最早由Slepian(1983)提出，旨在解決一維連續(xù)情況下的時域和頻域能量集中問題，隨后被不斷發(fā)展并引入到重力場研究領(lǐng)域中(Albertella et al., 1999; Simons and Dahlen, 2006; Simons et al., 2006; Simons, 2010).Slepian方法通過構(gòu)建在研究區(qū)域內(nèi)部的正交基函數(shù)(Simons, 2010; Slobbe et al., 2012; 朱廣彬等, 2012)，可使信號能量集中在研究區(qū)域內(nèi)部，適用于表達(dá)區(qū)域尺度的重力場變化，例如月球區(qū)域重力場模型(Han, 2008; 孫雪梅等, 2015)和中國大陸區(qū)域重力場模型(陳石等, 2017)，并且Slepian展開與相同階次的經(jīng)典球諧展開具有同樣的意義(Simons, 2010; Slobbe et al., 2012; 陳石等, 2017).

本文對利用Slepian局部譜分析恢復(fù)區(qū)域重力場的方法進(jìn)行了研究，并將Slepian方法應(yīng)用到我國華北地區(qū)，確定了華北年尺度的區(qū)域重力場時變模型.本文分以下3個部分展開研究，第1部分主要介紹Slepian方法的基本原理，第2部分首先說明觀測數(shù)據(jù)處理過程及重復(fù)測點(diǎn)分布，接著利用EGM2008模型和實(shí)際觀測點(diǎn)位分布，討論分析Slepian最佳截斷數(shù)和最大展開階數(shù)的確定，最后對多期地面實(shí)測重力觀測數(shù)據(jù)進(jìn)行處理，獲得我國華北地區(qū)2011—2013年時變重力場模型，并與研究區(qū)域內(nèi)GRACE衛(wèi)星重力估計(jì)結(jié)果進(jìn)行對比分析，第3部分給出討論與結(jié)論，并對下一步工作做出展望.

1 Slepian方法基本原理

(1)

(2)

其中，sn m為Slepian系數(shù).Slepian基函數(shù)在區(qū)域Γ內(nèi)的集中程度由(3)式給出

(3)

其中，κ為“聚集因子”(concentration factor)，取值0<κ<1.公式(3)的最大化可通過求解以下特征方程實(shí)現(xiàn)：

(4)

κ是上述特征方程的特征值，s是特征向量，即公式(2)中的Slepian系數(shù).矩陣D各元素由公式(5)給出

(5)

由公式(4)可解得(N+1)2個特征值和特征向量，將特征向量s代入式(2)，經(jīng)整理可得Slepian基函數(shù)sα[α=1, …, (N+1)2]具有整個球面Ω內(nèi)及區(qū)域Γ內(nèi)的正交特性：

(6)

其中δ為Kronecker 函數(shù).

(7)

其中，tα為Slepian展開系數(shù)，J為Slepian截斷數(shù).需要指出的是，對于球面區(qū)域Γ內(nèi)的重力場信號表達(dá)，Slepian展開與相同階次的經(jīng)典球諧展開具有同樣的意義(Simons, 2010; Slobbe et al., 2012; 陳石等, 2017).

公式(7)表明Slepian方法是一種具有稀疏(sparse)特性的時空局部譜方法.以本文研究的華北地區(qū)為例，圖1直觀反映了Slepian函數(shù)分別展開到90、120和180階時聚集因子κ的稀疏分布情況.對于120階Slepian展開，聚集因子κ的總數(shù)量為(120+1)×(120+1)=14641個，其中大于0.001的共41個，僅占總數(shù)量的0.3%，剩余14600個κ值均接近0，分布具有明顯的稀疏性.

圖1 華北地區(qū)對應(yīng)90、120和180階Slepian展開的聚集因子分布Fig.1 The concentration factors of Slepian functions in North China for N=90,120, and 180, respectively

圖2進(jìn)一步給出了華北地區(qū)Slepian展開到90、120和180階時基函數(shù)的空間分布情況，其中洋紅色曲線標(biāo)記重力數(shù)據(jù)覆蓋范圍的邊界，也即本文研究區(qū)域.由圖2可以看出，對于每組展開，當(dāng)聚集因子κ較接近1時，所對應(yīng)基函數(shù)基本都集中在研究區(qū)域內(nèi)部，能量很少泄漏到區(qū)域外部；當(dāng)聚集因子κ很小時(如圖2第3行中κ接近0.001的情況)，所對應(yīng)基函數(shù)基本集中在研究區(qū)域外部，在研究區(qū)域內(nèi)部幾乎沒有貢獻(xiàn).因此，在表達(dá)研究區(qū)域內(nèi)重力場信號時，可按公式(7)做截斷近似處理，只采用較大的聚集因子所對應(yīng)的基函數(shù).

圖2 華北地區(qū)3組Slepian展開的基函數(shù)空間分布，分別對應(yīng)90階(左邊2列)、120階(中間2列)和180階(右邊2列)展開.對于每組展開，第1和第2行為聚集因子κ最大的前4個基函數(shù)，第3行為κ接近0.001的基函數(shù)洋紅色線為研究區(qū)域邊界，黑色線為陸海邊界.Fig.2 Spatial patterns of the Slepian functions with corresponding concentration factors in North China for N=90 (the left 2 columns), N=120 (the central 2 columns), and N=180 (the right 2 columns). For each group, the first and second rows show the basis functions with the 4 largest concentration factors, the third row shows the 2 basis functions with concentration factors around the value of 0.001The magenta line denotes the boundary of the study area, while the black line indicates the land-ocean boundary.

2 數(shù)值計(jì)算和分析

2.1 研究區(qū)域及重復(fù)測點(diǎn)空間分布

本文采用了華北地區(qū)重力測網(wǎng)2011—2013年間的實(shí)際重力觀測成果，一共6期觀測數(shù)據(jù)(每年上下半年各測1期).重力網(wǎng)平差采用基于貝葉斯平差模型的pyBACGS軟件完成(Chen et al., 2019)，該方法可以較好地解決多臺相對重力儀器聯(lián)合平差的權(quán)系數(shù)優(yōu)化、非線性漂移估計(jì)和儀器格值的最優(yōu)化估計(jì)等問題，6期觀測平差結(jié)果的點(diǎn)值精度均優(yōu)于10×10-8m·s-2.圖3給出了6期數(shù)據(jù)覆蓋的重復(fù)測點(diǎn)位置的空間分布，共有364個.以最邊緣的重復(fù)觀測點(diǎn)位置確定閉合邊界(圖3中洋紅色曲線)，邊界內(nèi)部即為本文研究區(qū)域.由圖3可以看出，研究區(qū)域內(nèi)部重力測點(diǎn)的空間分布很不均勻：京津地區(qū)最為密集，該地區(qū)約占研究區(qū)域總面積的5.5%，測點(diǎn)數(shù)約占全部的15.9%，相對測點(diǎn)密度(測點(diǎn)百分比/面積百分比)約2.9；河南、安徽及山東境內(nèi)較為稀疏，該地區(qū)約占研究區(qū)域總面積的42.4%，測點(diǎn)數(shù)僅占全部的22.5%，相對測點(diǎn)密度約0.5.研究區(qū)域內(nèi)相鄰測點(diǎn)最小距離約0.2 km，最大距離約167.2 km，平均距離為43.5 km.

圖3 本文研究區(qū)域及重復(fù)測點(diǎn)位置分布洋紅色線為研究區(qū)域邊界，綠色十字代表重復(fù)測點(diǎn)位置，灰色線為省界，黑色線為陸海邊界.Fig.3 The study area and the locations of the repeated gravity stationsThe magenta line denotes the boundary of the study area, the green crosses show the locations of the repeated gravity stations, the gray lines show the provincial borders, and the black line indicates the land-ocean boundary.

在利用實(shí)測數(shù)據(jù)和Slepian方法恢復(fù)華北區(qū)域重力場前，還需解決兩個問題：(1) Slepian展開最佳截斷數(shù)的確定.研究區(qū)域內(nèi)觀測點(diǎn)僅有364個，應(yīng)用Slepian方法時，不管最大展開階數(shù)N如何取值，都須采用公式(7)給出的截斷形式(J≤364)，否則求解Slepian展開系數(shù)時會出現(xiàn)欠定情況；(2) Slepian最大展開階數(shù)N的確定.研究區(qū)域內(nèi)觀測點(diǎn)分布稀疏不均，需測試基于現(xiàn)有測點(diǎn)分布對重力變化的恢復(fù)能力，以此確定Slepian的最大展開階數(shù)N.接下來將在2.2和2.3節(jié)對以上兩個問題進(jìn)行詳細(xì)討論.

2.2 Slepian展開最佳截斷數(shù)的確定

對于不同的Slepian展開階數(shù)，由于聚集因子和基函數(shù)的分布不同(圖1和2)，在根據(jù)公式(7)確定重力場模型時選取的截斷數(shù)也有所不同.假設(shè)觀測點(diǎn)數(shù)為M，最佳截斷數(shù)Jopt會在Msha和M之間變化，Msha為Shannon數(shù)，即最小截斷數(shù)(Simons, 2010; Slobbe et al., 2012)

(8)

其中，κα為第α個聚集因子，N為最大展開階數(shù)，A/4π為研究區(qū)域所占球面面積的比例.

對Slepian最佳截斷數(shù)Jopt的確定，目前比較常用的方法為：首先給定一個已知N階的初始場，然后模擬測點(diǎn)處的觀測值，接著利用模擬觀測值和Slepian方法基于變化的J值(Msha≤J≤M)獲得一系列N階恢復(fù)結(jié)果，最后將恢復(fù)結(jié)果與初始場進(jìn)行比較，令二者差值均方根(root mean square, RMS)最小的J即為Slepian展開最佳截斷數(shù)Jopt(Slobbe et al., 2012; Plattner and Simons, 2017; 陳石等, 2017).這一確定準(zhǔn)則雖然簡便易行，但是(1)采用的已知場跟恢復(fù)結(jié)果的階數(shù)是一致的，不適合二者階數(shù)不同的情形.本文的解算數(shù)據(jù)為地面重力觀測值，信號為全波段(可展開至無窮階)，遠(yuǎn)大于所要恢復(fù)的N階；(2)沒有顧及模型的復(fù)雜度，在實(shí)際應(yīng)用中可能存在模型參數(shù)過多導(dǎo)致的過擬合現(xiàn)象.

對于確定Jopt時存在的問題(1)，應(yīng)用Slepian理論求解N階模型時，觀測數(shù)據(jù)中高于N階的信號會影響解算得到的低階系數(shù)，造成寬帶泄漏(broadband leakage)(Simons and Dahlen, 2006; Slobbe et al., 2012).為削弱寬帶泄漏的影響，需盡可能地移除觀測數(shù)據(jù)中高于N階的信號.本文采用Bomfim等(2013)在比較地面重力觀測數(shù)據(jù)和GOCE(Gravity Field and Steady-State Ocean Circulation Explorer)衛(wèi)星重力數(shù)據(jù)時采用的濾波方法，移除地面重力數(shù)據(jù)中高于N階的高頻信號.

對于問題(2)，本文采用貝葉斯信息量準(zhǔn)則(Bayesian Information Criterion, BIC)(Baur, 2012; He et al., 2019)確定最佳截斷數(shù)Jopt.通過引入關(guān)于參數(shù)數(shù)量的懲罰項(xiàng)，BIC在提高擬合度的同時，要求盡可能少的模型參數(shù)以避免過擬合現(xiàn)象(Baur, 2012; He et al., 2019).對于截斷數(shù)為J的恢復(fù)結(jié)果，

BIC=J·ln(M)-2ln(L)，

(9)

其中，M為觀測值個數(shù)，L為最大似然函數(shù)值.當(dāng)BIC值最小時，對應(yīng)的J即為最佳截斷數(shù)Jopt，此時應(yīng)同樣滿足Msha≤Jopt≤M.

本文利用BIC準(zhǔn)則確定了華北地區(qū)90、120和180階(半波長空間分辨率約200 km、150 km和100 km)Slepian展開時的最佳截斷數(shù)，具體步驟如下：首先選取EGM2008重力場模型(Pavlis et al., 2012)作為已知場，以模型最高2190階自由空氣重力異常模擬全波段的地面重力觀測值，獲得研究區(qū)域內(nèi)地面重復(fù)測點(diǎn)處重力值.然后利用模擬觀測值和Slepian方法恢復(fù)N階重力場，最后將對應(yīng)不同J值(2≤J≤364)的恢復(fù)結(jié)果與初始N階EGM2008模型進(jìn)行比較，基于公式(9)確定對應(yīng)的BIC值.上述過程執(zhí)行3次，N分別取90、120和180，獲得3組BIC值隨J值的變化趨勢(圖4).從J=2開始，每組BIC曲線首先下降到達(dá)最小值，然后呈上升趨勢；當(dāng)J>50時，3組BIC曲線均為上升趨勢，因此圖4僅列出了2≤J≤121的部分.對于華北地區(qū)90、120和180階Slepian展開，最大截斷數(shù)為M=364(即總測點(diǎn)數(shù))，最小截斷數(shù)(Shannon數(shù))Msha分別為8、14和32，BIC曲線的最小值分別出現(xiàn)在J=12、J=26和J=36處(圖4)，即以上90、120和180階Slepian展開的最佳截斷數(shù)Jopt分別為12、26和36.

圖4 華北地區(qū)BIC與J取值的關(guān)系圖(N=90,120,180)，虛線指示BIC最小值對應(yīng)的J值Fig.4 The relationship between BIC and J in North China (N=90,120,180), where a dashed line indicates the smallest BIC and the corresponding J value in each scenario

表1給出了N取90、120和180時，分別基于Msha和BIC確定的Jopt進(jìn)行截斷，Slepian方法對華北地區(qū)初始模型(即N階EGM2008模型)恢復(fù)情況的統(tǒng)計(jì)信息.

表1 選取不同截斷數(shù)時Slepian恢復(fù)結(jié)果在研究區(qū)域內(nèi)評估信息Table 1 Quality assessment of the Slepian-recovered results based on different truncation numbers within the study area

圖5以120階展開為例，給出了基于不同截斷數(shù)恢復(fù)結(jié)果的空間分布及其與初始模型的差異(格網(wǎng)間隔為0.1°×0.1°).結(jié)合表1和圖5可以發(fā)現(xiàn)，采用BIC確定的Jopt截斷較使用Msha截斷，Slepian恢復(fù)結(jié)果有較為明顯的改進(jìn)：(1)基于不同N的3組測試，以BIC確定的Jopt截斷時，Slepian恢復(fù)結(jié)果的精度均有較顯著提升；(2)基于Jopt和Msha的120階恢復(fù)結(jié)果在空間分布上較為接近，但與初始EGM2008模型做差后，基于Jopt恢復(fù)結(jié)果的殘差在山東省西部和安徽省西北部明顯減小，表明以Jopt進(jìn)行截斷的恢復(fù)結(jié)果更加接近初始模型.基于90階和180階恢復(fù)結(jié)果的空間分布可得出同樣結(jié)論，相關(guān)結(jié)果不再逐一列出.

圖5 華北地區(qū)120階Slepian恢復(fù)結(jié)果及其與初始EGM2008模型的差異(a) 以Shannon數(shù)Msha進(jìn)行截斷的恢復(fù)結(jié)果; (b) 以BIC準(zhǔn)則確定的最佳截斷數(shù)Jopt進(jìn)行截斷的恢復(fù)結(jié)果; (c) 基于Msha的恢復(fù)結(jié)果與EGM2008模型的差值; (d) 基于Jopt的恢復(fù)結(jié)果與EGM2008模型的差值.洋紅色線為研究區(qū)域邊界，灰色線為省界，黑色線為陸海邊界.Fig.5 The Slepian recovered results (degree 120) and differences with the EGM2008 inputs(a) Slepian results after truncating at Msha; (b) Slepian results after truncating at BIC-derived Jopt; (c) Recovered results in (a) minus the EGM2008 inputs; (d) Recovered results in (b) minus the EGM2008 inputs. The magenta line denotes the boundary of the study area, the gray lines show the provincial borders, and the black line indicates the land-ocean boundary.

2.3 Slepian最大展開階數(shù)的確定

本節(jié)將進(jìn)一步討論基于現(xiàn)有測點(diǎn)分布和本文方法所能恢復(fù)研究區(qū)域內(nèi)重力場變化的最大階數(shù).根據(jù)本文2.1節(jié)中統(tǒng)計(jì)信息，測點(diǎn)稀疏區(qū)域內(nèi)相鄰測點(diǎn)的最大距離約167 km，較為接近120階重力場模型的半波長空間分辨率(約150 km)，明顯小于180階的空間分辨率(約100 km)，但僅從測點(diǎn)最大距離難以確定最大展開階數(shù).本文2.2節(jié)基于研究區(qū)域內(nèi)現(xiàn)有測點(diǎn)分布和Slepian方法進(jìn)行重復(fù)實(shí)驗(yàn)，給出了Slepian方法對研究區(qū)域內(nèi)90、120和180階EGM2008初始模型恢復(fù)情況的對比分析(表1)，結(jié)果表明120階恢復(fù)結(jié)果與EGM2008初始場差值的RMS最小，即120階恢復(fù)結(jié)果的精度最高.基于上述對比分析，本文選取研究區(qū)域內(nèi)Slepian最大展開階數(shù)為N=120.

2.4 觀測點(diǎn)位分布對Slepian恢復(fù)結(jié)果的影響

基于實(shí)際重力觀測恢復(fù)區(qū)域重力場時，Slepian恢復(fù)結(jié)果受到寬頻泄漏和地表重力觀測點(diǎn)分布不均勻(疏密)程度的影響，例如2.2節(jié)中的實(shí)驗(yàn)結(jié)果(表1和圖5).本節(jié)將著重討論華北地區(qū)現(xiàn)有測點(diǎn)分布對恢復(fù)研究區(qū)域內(nèi)部120階重力場變化的影響.選取研究區(qū)域內(nèi)部以0.1°×0.1°間隔均勻分布的5019個點(diǎn)作為“控制點(diǎn)”，通過比較“控制點(diǎn)”處的EGM2008模型值與Slepian恢復(fù)值，對Slepian方法解算的120階模型進(jìn)行評估.為了避免寬頻泄漏的影響，采用120階EGM2008重力場模型模擬地面觀測值(此時Jopt=21由BIC準(zhǔn)則確定).表2列出了研究區(qū)域內(nèi)相關(guān)評估結(jié)果，在研究區(qū)域內(nèi)EGM2008已知場重力變化約58×10-5m·s-2，Slepian解算模型的變化約55×10-5m·s-2，與已知場吻合較好，整體偏差較小(-0.13×10-5m·s-2)，解算誤差為2.50×10-5m·s-2.模型解算誤差受地表重力觀測點(diǎn)分布不均勻(疏密)程度的影響較為顯著，在觀測點(diǎn)位最為密集的京津地區(qū)(“控制點(diǎn)”數(shù)為287個)，相對測點(diǎn)密度約2.9(相對測點(diǎn)定義見2.1節(jié))，解算模型在此區(qū)域的誤差明顯減小，為2.05×10-5m·s-2;在觀測點(diǎn)位較為稀疏的河南、安徽和山東(“控制點(diǎn)”數(shù)為2071個)，相對測點(diǎn)密度約0.5，解算模型在此區(qū)域的誤差顯著增大至2.98×10-5m·s-2.未來若能在河南、安徽和山東加密重力測網(wǎng)、增加更多的重復(fù)重力觀測點(diǎn)，將較大地提高華北區(qū)域重力場的解算精度.

表2 觀測點(diǎn)位分布對Slepian恢復(fù)結(jié)果(N=120)的影響(單位：10-5 m·s-2)Table 2 Impact of the distribution of gravity stations on the Slepian-recovered results (N=120)(unit: 10-5 m·s-2)

2.5 基于地面實(shí)測數(shù)據(jù)的華北2011—2013時變重力場

在研究區(qū)域內(nèi)部364個重復(fù)測點(diǎn)處，將相鄰兩年對應(yīng)測期內(nèi)的重力觀測做差(后一年觀測減去前一年)獲得1年尺度重力變化.所有比較都限定在相對應(yīng)的測期內(nèi)，以減弱水文變化等季節(jié)性信號的影響(祝意青等,2012,2013; 徐偉民等,2016).最終得到4組1年尺度重力變化，即基于第1期觀測的2013-04—2012-03、2012-03—2011-03(“—”號表示前后兩組數(shù)據(jù)做差)，以及基于第2期觀測的2013-08—2012-08和2012-08—2011-08.獲得研究區(qū)域重復(fù)測點(diǎn)處的重力變化后，接下來應(yīng)用Slepian方法恢復(fù)華北地區(qū)120階時變重力場模型，空間變化結(jié)果(格網(wǎng)間距為0.1°×0.1°)如圖6所示.

圖6反映了2011—2013年間華北地區(qū)較為明顯的重力變化信號.圖6a中大同—太原一帶正重力變化，圖6c中石家莊—北京一帶的正重力變化、太原—鄭州—合肥一帶的負(fù)重力變化以及圖6d中大范圍的正重力變化，與徐偉民等(2016)給出的結(jié)果基本一致.除Slepian解算結(jié)果外，本文還利用最小曲率插值法(Wessel et al., 2019)得到了2013-04—2012-03、2012-03—2011-03、2013-08—2012-08和2012-08—2011-08等4組重力變化的空間插值結(jié)果(格網(wǎng)間距為0.1°×0.1°)，如圖7所示.

圖6 華北地區(qū)由Slepian方法確定的120階1年尺度重力場變化(2011—2013)空間分布(a) (c) 基于相鄰兩年第1期觀測差值; (b) (d) 基于相鄰兩年第2期觀測差值.洋紅色線為研究區(qū)域邊界，灰色線為省界，黑色線為陸海邊界.Fig.6 Spatial pattern of the annual gravity changes up to degree 120 determined by the Slepian method in North China from 2011 to 2013(a) (c) Based on the differences between observations from the first campaigns of two adjacent years; (b) (d) Based on the differences between observations from the second campaigns of two adjacent years. The magenta line denotes the boundary of the study area, the gray lines show the provincial borders, and the black line indicates the land-ocean boundary.

圖7 華北地區(qū)1年尺度重力場變化(2011—2013)最小曲率空間插值結(jié)果(a) (c) 基于相鄰兩年第1期觀測差值; (b) (d) 基于相鄰兩年第2期觀測差值.洋紅色線為研究區(qū)域邊界，灰色線為省界，黑色線為陸海邊界.Fig.7 Spatial pattern of the annual gravity changes obtained by minimum curvature interpolation in North China from 2011 to 2013(a) (c) Based on the differences between observations from the first campaigns of two adjacent years; (b) (d) Based on the differences between observations from the second campaigns of two adjacent years. The magenta line denotes the boundary of the study area, the gray lines show the provincial borders, and the black line indicates the land-ocean boundary.

對比圖6和圖7可以發(fā)現(xiàn)，Slepian方法所確定的模型重力信號基本都集中在研究區(qū)域內(nèi)部且較為光滑，可有效壓制高頻干擾，空間插值方法所得結(jié)果存在較多的高頻信號.值得強(qiáng)調(diào)的是，Slepian方法在頻率域內(nèi)恢復(fù)重力場，可通過調(diào)整展開階數(shù)從觀測數(shù)據(jù)中提取所需的信號頻段，這一點(diǎn)利用空間插值方法無法實(shí)現(xiàn).

除1年尺度時變重力場模型外，本文還解算了2組華北地區(qū)120階的2年尺度重力變化模型(最優(yōu)截斷數(shù)Jopt同樣取26)，第1組利用基于上半年觀測(2013-04—2011-03)的重力變化數(shù)據(jù)解算，第2組利用基于下半年觀測(2013-08—2011-08)的重力變化數(shù)據(jù)解算，解算結(jié)果如圖8所示.

圖8 華北地區(qū)由Slepian方法確定的120階2年尺度重力場變化(2011—2013)空間分布(a) 基于2013-04和2011-03的觀測差值; (b) 基于2013-08和2011-08的觀測差值. 洋紅色線為研究區(qū)域邊界，灰色線為省界，黑色線為陸海邊界.Fig.8 Spatial pattern of the 2-year gravity changes (between 2011 and 2013) up to degree 120 determined by the Slepian method in North China(a) From differences between 2013-04 and 2011-03; (b) From differences between 2013-08 and 2011-08. The magenta line denotes the boundary of the study area, the gray lines show the provincial borders, and the black line indicates the land-ocean boundary.

在研究區(qū)域內(nèi)部，圖8a和圖8b所示兩組120階結(jié)果在變化極值上存在較大差異，圖8a結(jié)果最大值和最小值分別為16.94×10-8m·s-2、-79.04×10-8m·s-2，而圖8b結(jié)果最大值和最小值分別為52.68×10-8m·s-2、-41.70×10-8m·s-2.除去重力變化極值差異外，圖8a和圖8b所示重力變化在空間分布上大致呈現(xiàn)“南減北增”的特征：在河南省南部和安徽省西北部(鄭州—合肥連線以南區(qū)域)，重力變化呈減少趨勢；在河北省南部、京津和大同附近區(qū)域，重力變化呈增加趨勢.

2.6 基于GRACE數(shù)據(jù)的華北2011—2013時變重力場

為了與2.5節(jié)中地面結(jié)果進(jìn)行對比，本文利用Slepian方法和GRACE 衛(wèi)星數(shù)據(jù)估計(jì)了華北地區(qū)2011—2013年間的1年尺度重力變化，結(jié)果如圖9所示.本文采用由美國德克薩斯大學(xué)空間研究中心(Center for Space Research at the University of Texas)發(fā)布的GRACE Level-2 Release 06 (RL06)月時變重力場模型，數(shù)據(jù)為球諧模型形式，完整展開到60階.由于GRACE無法觀測地心運(yùn)動(球諧系數(shù)中相關(guān)1階項(xiàng)數(shù)值均為0)，同時2階項(xiàng)系數(shù)C20也包含較大誤差，本文對RL06球諧系數(shù)做了如下替換：將原始1階項(xiàng)C10、C11和S11替換為Sun等(2016)基于大氣海洋模型、冰后回彈模型和GRACE數(shù)據(jù)計(jì)算的結(jié)果，將原始C20項(xiàng)替換為衛(wèi)星激光測距確定的結(jié)果(Cheng and Ries, 2017).

圖9 華北地區(qū)由Slepian方法得到的1年尺度GRACE重力變化(2011—2013)空間分布(a),(b),(c),(d)分別為2012-03—2011-03、2012-08—2011-08、2013-04—2012-03和2013-07—2012-08結(jié)果. 洋紅色線為研究區(qū)域邊界，灰色線為省界，黑色線為陸海邊界.Fig.9 Spatial pattern of the annual gravity changes from GRACE obtained by the Slepian method in North China (2011—2013)(a), (b), (c) and (d) denote results of 2012-03—2011-03, 2012-08—2011-08, 2013-04—2012-03, and 2013-07—2012-08, respectively. The magenta line denotes the boundary of the study area, the gray lines show the provincial borders, and the black line indicates the land-ocean boundary.

Slepian方法由Harig和Simons(2012)引入到GRACE數(shù)據(jù)后處理中，首先將GRACE球諧系數(shù)轉(zhuǎn)換為Slepian系數(shù)，然后利用Slepian系數(shù)計(jì)算重力場變化(Harig and Simons, 2012, 2015; Von Hippel and Harig, 2019; 高春春等, 2019).本文在RL06球諧系數(shù)的轉(zhuǎn)換過程中，根據(jù)2.2節(jié)方法將60階Slepian展開的最佳截斷數(shù)確定為Jopt=5(Shannon數(shù)為4).為了跟圖6中地面時變重力觀測時段保持一致， 本文盡量選取GRACE相同時段數(shù)據(jù)；對于地面重力的2013-08—2012-08時間段，由于GRACE缺失2013-08和2013-09月數(shù)據(jù)，因此替換為最接近的2013-07—2012-08.利用Slepian方法完成對GRACE數(shù)據(jù)的處理后，本文未再扣除地殼均衡調(diào)整(Glacial Isostatic Adjustment)和水文變化的影響.

本文還利用Slepian方法計(jì)算了華北地區(qū)2011—2013年間2年尺度的GRACE重力變化(2013-04—2011-03和2013-07—2011-08)，如圖10所示.

圖10 華北地區(qū)由Slepian方法得到的2年尺度GRACE重力變化空間分布(a)，(b)分別為2013-04—2011-03和2013-07—2011-08結(jié)果. 洋紅色線為研究區(qū)域邊界，灰色線為省界，黑色線為陸海邊界.Fig.10 Spatial pattern of the 2-year gravity changes from GRACE derived by the Slepian method in North China(a) From differences between 2013-04 and 2011-03; (b) From differences between 2013-07 and 2011-08. The magenta line denotes the boundary of the study area, the gray lines show the provincial borders, and the black line indicates the land-ocean boundary.

2.7 地面和GRACE結(jié)果的比較

利用Slepian方法獲得華北地區(qū)2011—2013基于地面和GRACE觀測的重力場變化后，本節(jié)就1年和2年尺度重力變化的空間分布對二者進(jìn)行了比較.總體上看，基于地面重力觀測的120階結(jié)果(空間分辨率約150 km)給出的重力變化量級更大、細(xì)節(jié)(高頻信號)更多.由于空間分辨率相差較大，基于GRACE衛(wèi)星觀測的結(jié)果(60階，空間分辨率約400 km)無法給出與地面120階結(jié)果相當(dāng)?shù)闹亓ψ兓畔ⅲ咴跁r空變化的顯著趨勢上存在較好的對應(yīng).

對于1年尺度重力變化，首先對比時間段不完全一致的2013-08—2012-08地面結(jié)果和2013-07—2012-08 GRACE結(jié)果，前者重力變化在研究區(qū)域內(nèi)呈整體增加趨勢(圖6d)，后者大致呈“北減南增”趨勢(圖9d)，二者相差較大.當(dāng)?shù)孛娼Y(jié)果與GRACE結(jié)果所取時間段完全一致時，二者出現(xiàn)顯著變化的區(qū)域大致相同，如2012-03—2011-03(圖6a和圖9a)、2012-08—2011-08(圖6b和圖9b)和2013-04—2012-03(圖6c和圖9c).以2012-08—2011-08時間段為例，地面結(jié)果在河南與安徽西北部呈明顯減少趨勢(圖6b)，GRACE結(jié)果在河南和安徽也呈明顯減少趨勢，由于分辨率原因，GRACE重力顯著減小區(qū)域的范圍更大、中心也偏北一些(圖9b).

對于2年尺度的重力變化，地面結(jié)果(圖8)與GRACE結(jié)果(圖10)表現(xiàn)出較為一致的“南減北增”空間分布特征：2013-04—2011-03時間段，GRACE結(jié)果在研究區(qū)域內(nèi)部基本為減少趨勢，在京津地區(qū)呈增加趨勢(圖10a)，地面結(jié)果在山西省南部、河南省、安徽省西北部及山東省西南部呈重力減少趨勢，在京津地區(qū)、河北省以及山西省北部(大同鄰近區(qū)域)呈增加趨勢(圖8a)；2013-08—2011-08(衛(wèi)星數(shù)據(jù)為2013-07—2011-08)時間段，GRACE結(jié)果大致以太原—石家莊—衡水連線為分界，以南區(qū)域呈減少趨勢、以北區(qū)域呈增加趨勢(圖10b)，地面結(jié)果重力減少區(qū)域出現(xiàn)在鄭州—合肥連線以南，以北區(qū)域呈差異性增加趨勢，其中山東省西南部增加幅度較小，山西省南部稍大，京津地區(qū)、河北省以及山西省北部最大(圖8b).

進(jìn)一步對比圖8和圖10，可發(fā)現(xiàn)2年尺度的地面與GRACE衛(wèi)星結(jié)果差異最大區(qū)域出現(xiàn)在石家莊—衡水連線一帶，2011—2013年區(qū)域內(nèi)地面結(jié)果為增大趨勢，GRACE結(jié)果為減少趨勢.該區(qū)域長期存在嚴(yán)重的淺層、深層地下水超量開采以及地表沉降問題(張永紅等, 2016; 馮偉等, 2017; Zhao et al., 2019)，這是影響區(qū)域內(nèi)重力場長期變化的兩個主要因素(Shen et al., 2015).對地面測量而言，地下水虧損會令重力觀測減小，地表沉降會導(dǎo)致重力觀測增加.根據(jù)Shen等(2015)的地面觀測結(jié)果，2009—2013年在河北省南部(包含石家莊—衡水連線一帶)由地表沉降導(dǎo)致的重力增加效應(yīng)要大于地下水虧損帶來的減小效應(yīng)；由于2009—2013年間石家莊—衡水連線一帶地表沉降速度基本維持不變(約-0.05 m·a-1, Zhao et al., 2019)，我們推斷2011—2013年區(qū)域內(nèi)地表沉降的影響仍然大于地下水虧損的影響，因此地面重力結(jié)果呈增大趨勢.對GRACE觀測而言，在衛(wèi)星軌道處地面沉降導(dǎo)致的重力變化已大幅衰減，此時地下水虧損導(dǎo)致的重力減小效應(yīng)占主導(dǎo)地位，因此GRACE結(jié)果整體為減少趨勢，馮偉等(2017)給出的2011、2012和2013年GRACE結(jié)果在該區(qū)域內(nèi)的減少趨勢也證實(shí)了這一點(diǎn).

3 討論與結(jié)論

Slepian方法的基函數(shù)除具備整個球面上的正交性之外，同時具備在指定區(qū)域內(nèi)部的正交性，適用于構(gòu)建局部重力場模型.通過構(gòu)建在研究區(qū)域內(nèi)部的正交基函數(shù)，可使信號能量集中在區(qū)域內(nèi)部，從而有效減小信號泄漏的影響，提高研究區(qū)域內(nèi)恢復(fù)結(jié)果的信噪比(Harig and Simons, 2012, 2015; 高春春等, 2019)，這一特點(diǎn)利于提取和分析較小區(qū)域內(nèi)的重力信號.

Slepian基函數(shù)另一個重要的特點(diǎn)即具備稀疏性(Simons, 2010)：以120階Slepian展開為例，基函數(shù)共有(120+1)2=14641個，其中聚集因子大于0.001的僅有41個(占總數(shù)量的0.3%)，而這41個基函數(shù)的累積能量占全部Slepian基函數(shù)能量的99.98%.根據(jù)Slepian基函數(shù)具備稀疏性這一特點(diǎn)，可對Slepian表達(dá)進(jìn)行截斷、只采用聚集因子較大的若干項(xiàng)基函數(shù)來表達(dá)研究區(qū)域內(nèi)的重力場信號.本文2.2節(jié)基于BIC準(zhǔn)則確定最佳截斷數(shù)Jopt對Slepian展開進(jìn)行近似(即只采用前Jopt項(xiàng)基函數(shù))，以ε1表示近似誤差(包含寬頻泄漏誤差)，當(dāng)展開階數(shù)取120時，ε1約為初始信號大小的9%.由于Slepian展開與相同階次的經(jīng)典球諧展開具有同樣的意義(Simons, 2010; Slobbe et al., 2012; 陳石等, 2017)，當(dāng)將重力觀測數(shù)據(jù)展開到N階且Slepian方法取其全部(N+1)2個基函數(shù)時，二者的“截斷誤差”ε2(重力信號在N+1～∞階之間的部分)是一致的；若Slepian方法取Jopt個基函數(shù)時，Slepian展開的“截斷誤差”除ε2外，還要包含近似誤差ε1(Han and Razeghi, 2017).

稀疏性使得利用Slepian方法恢復(fù)區(qū)域重力場時所需的重力觀測值數(shù)量大為減少，仍以構(gòu)建120階重力場為例，在不考慮測點(diǎn)分布的情況下，至少需要14641個觀測值才能求解全部展開系數(shù)；而根據(jù)稀疏性，需要最少Jopt個觀測值(Jopt遠(yuǎn)小于 14641)即可求解所需的展開系數(shù).

本文研究了基于Slepian譜分析理論的局部重力場建模方法，并利用EGM2008模型和實(shí)際觀測點(diǎn)位分布，討論分析了Slepian最佳截斷數(shù)和最大展開階數(shù)的確定，最后利用Slepian方法和2011—2013年間地面實(shí)測重力數(shù)據(jù)對華北區(qū)域時變重力場構(gòu)建做出了嘗試，獲得了華北地區(qū)120階(約150 km空間分辨率)的2011—2013年變化重力場模型，并與研究區(qū)域內(nèi)基于Slepian方法和GRACE衛(wèi)星數(shù)據(jù)得到的重力變化進(jìn)行了比較.研究結(jié)果表明：

(1)貝葉斯信息量準(zhǔn)則(BIC)可作為確定Slepian展開最佳截斷數(shù)的有效手段.對于90、120和180階Slepian展開(半波長空間分辨率約200 km、150 km和100 km)，基函數(shù)數(shù)量分別為8281、14641和32761個，基于BIC的最佳截斷數(shù)分別為12、26和36，表明Slepian方法是一種具有稀疏特性的時空局部譜方法.

(2)對于本文選取的3組展開階數(shù)(90、120和180)，研究區(qū)域內(nèi)Slepian展開的最大階數(shù)為120，模擬實(shí)驗(yàn)的解算誤差為4.98×10-5m·s-2(包含測點(diǎn)分布和寬頻泄漏的影響).僅考慮測點(diǎn)分布對Slepian展開的影響時，120階模擬實(shí)驗(yàn)的解算誤差變?yōu)?.50×10-5m·s-2; 在重復(fù)測點(diǎn)最為密集的京津地區(qū)，模型解算誤差為2.05×10-5m·s-2; 在重復(fù)測點(diǎn)稀疏的河南、安徽和山東區(qū)域，解算誤差增大至2.98×10-5m·s-2.若能在重復(fù)測點(diǎn)稀疏地區(qū)加密重力測網(wǎng)，可較大提高華北區(qū)域重力場模型的解算精度.

(3)基于實(shí)測數(shù)據(jù)和Slepian方法解算的120階重力場模型，重力變化信號基本都集中在研究區(qū)域內(nèi)部，信號向研究區(qū)域外部的泄漏較??；可提供準(zhǔn)確的信號頻段，并濾除高頻干擾，凸顯主要的重力變化信號.

(4)2011—2013年，研究區(qū)域內(nèi)GRACE估計(jì)結(jié)果(60階)與地面結(jié)果(120階)在時空分布的顯著趨勢上存在較好的對應(yīng)，證明了本文利用Slepian方法和地面數(shù)據(jù)所得重力場模型的可靠性.此外，120階地面結(jié)果的空間差異性更強(qiáng)、重力變化的細(xì)節(jié)(高頻信號)更多，可作為GRACE衛(wèi)星模型的有效補(bǔ)充.

本文確定的華北地區(qū)120階2011—2013時變重力場模型數(shù)據(jù)可由https:∥github.com/jchhan/TVGM-NorthChina下載.基于本文的理論方法框架和研究成果，未來如果使用華北區(qū)域包含更多測期和測點(diǎn)的地面數(shù)據(jù)，可在時間跨度和空間分辨率上實(shí)現(xiàn)對本文2011—2013年變化結(jié)果的改進(jìn)，為華北構(gòu)造活動分析、地震風(fēng)險性評估以及水資源變化監(jiān)測等研究提供持續(xù)可靠的時變重力場模型產(chǎn)品.此外，基于本文思路還可以實(shí)現(xiàn)在Slepian域內(nèi)融合來自地面和衛(wèi)星觀測在內(nèi)的多源重力數(shù)據(jù)，有望獲得精度更高、適用范圍更廣的區(qū)域重力場模型.

致謝感謝三位評審專家提出的寶貴修改建議.感謝美國普林斯頓大學(xué)Frederik Simons教授提供Slepian局部譜分析程序包(https:∥github.com/fjsimons)和德國地學(xué)中心(GeoForschungsZentrum)提供GRACE RL06 Level-2數(shù)據(jù)(ftp:∥isdcftp.gfz-potsdam.de/grace/).