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

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

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

0 引言

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

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

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

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

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

1 Slepian方法基本原理

(1)

(2)

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

(3)

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

(4)

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

(5)

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

(6)

其中δ為Kronecker 函數.

(7)

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

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

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

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

圖2 華北地區3組Slepian展開的基函數空間分布，分別對應90階(左邊2列)、120階(中間2列)和180階(右邊2列)展開.對于每組展開，第1和第2行為聚集因子κ最大的前4個基函數，第3行為κ接近0.001的基函數洋紅色線為研究區域邊界，黑色線為陸海邊界.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 數值計算和分析

2.1 研究區域及重復測點空間分布

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

圖3 本文研究區域及重復測點位置分布洋紅色線為研究區域邊界，綠色十字代表重復測點位置，灰色線為省界，黑色線為陸海邊界.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.

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

2.2 Slepian展開最佳截斷數的確定

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

(8)

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

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

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

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

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

(9)

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

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

圖4 華北地區BIC與J取值的關系圖(N=90,120,180)，虛線指示BIC最小值對應的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進行截斷，Slepian方法對華北地區初始模型(即N階EGM2008模型)恢復情況的統計信息.

表1 選取不同截斷數時Slepian恢復結果在研究區域內評估信息Table 1 Quality assessment of the Slepian-recovered results based on different truncation numbers within the study area

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

圖5 華北地區120階Slepian恢復結果及其與初始EGM2008模型的差異(a) 以Shannon數Msha進行截斷的恢復結果; (b) 以BIC準則確定的最佳截斷數Jopt進行截斷的恢復結果; (c) 基于Msha的恢復結果與EGM2008模型的差值; (d) 基于Jopt的恢復結果與EGM2008模型的差值.洋紅色線為研究區域邊界，灰色線為省界，黑色線為陸海邊界.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最大展開階數的確定

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

2.4 觀測點位分布對Slepian恢復結果的影響

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

表2 觀測點位分布對Slepian恢復結果(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 基于地面實測數據的華北2011—2013時變重力場

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

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

圖6 華北地區由Slepian方法確定的120階1年尺度重力場變化(2011—2013)空間分布(a) (c) 基于相鄰兩年第1期觀測差值; (b) (d) 基于相鄰兩年第2期觀測差值.洋紅色線為研究區域邊界，灰色線為省界，黑色線為陸海邊界.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 華北地區1年尺度重力場變化(2011—2013)最小曲率空間插值結果(a) (c) 基于相鄰兩年第1期觀測差值; (b) (d) 基于相鄰兩年第2期觀測差值.洋紅色線為研究區域邊界，灰色線為省界，黑色線為陸海邊界.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可以發現，Slepian方法所確定的模型重力信號基本都集中在研究區域內部且較為光滑，可有效壓制高頻干擾，空間插值方法所得結果存在較多的高頻信號.值得強調的是，Slepian方法在頻率域內恢復重力場，可通過調整展開階數從觀測數據中提取所需的信號頻段，這一點利用空間插值方法無法實現.

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

圖8 華北地區由Slepian方法確定的120階2年尺度重力場變化(2011—2013)空間分布(a) 基于2013-04和2011-03的觀測差值; (b) 基于2013-08和2011-08的觀測差值. 洋紅色線為研究區域邊界，灰色線為省界，黑色線為陸海邊界.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.

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

2.6 基于GRACE數據的華北2011—2013時變重力場

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

圖9 華北地區由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結果. 洋紅色線為研究區域邊界，灰色線為省界，黑色線為陸海邊界.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數據后處理中，首先將GRACE球諧系數轉換為Slepian系數，然后利用Slepian系數計算重力場變化(Harig and Simons, 2012, 2015; Von Hippel and Harig, 2019; 高春春等, 2019).本文在RL06球諧系數的轉換過程中，根據2.2節方法將60階Slepian展開的最佳截斷數確定為Jopt=5(Shannon數為4).為了跟圖6中地面時變重力觀測時段保持一致， 本文盡量選取GRACE相同時段數據；對于地面重力的2013-08—2012-08時間段，由于GRACE缺失2013-08和2013-09月數據，因此替換為最接近的2013-07—2012-08.利用Slepian方法完成對GRACE數據的處理后，本文未再扣除地殼均衡調整(Glacial Isostatic Adjustment)和水文變化的影響.

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

圖10 華北地區由Slepian方法得到的2年尺度GRACE重力變化空間分布(a)，(b)分別為2013-04—2011-03和2013-07—2011-08結果. 洋紅色線為研究區域邊界，灰色線為省界，黑色線為陸海邊界.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結果的比較

利用Slepian方法獲得華北地區2011—2013基于地面和GRACE觀測的重力場變化后，本節就1年和2年尺度重力變化的空間分布對二者進行了比較.總體上看，基于地面重力觀測的120階結果(空間分辨率約150 km)給出的重力變化量級更大、細節(高頻信號)更多.由于空間分辨率相差較大，基于GRACE衛星觀測的結果(60階，空間分辨率約400 km)無法給出與地面120階結果相當的重力變化信息，但二者在時空變化的顯著趨勢上存在較好的對應.

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

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

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

3 討論與結論

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

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

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

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

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

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

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

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

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

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