利用改進的PnMl方法處理GRACE南北條帶誤差

1 數據和處理方法

1.1 數據

本研究采用3大機構(CSR、JPL、GFZ)2002-04～2019-09的GRACE/GFO 月時變重力場模型(即球諧系數模型)估計全球質量場變化。將模型截斷至60階次并移除背景重力場(平均重力場),獲得月重力場變化量。球諧系數模型的一階項由GRACE/GFO數據結合海洋數值模型確定[5],C20和C30項由相應的高精度衛星激光測距結果替代[6-7],冰川均衡改正則采用ICE-6G_D模型[8]。此外,采用CSR發布的RL06版本GRACE/GFO 0.25°×0.25° Mascon產品(CSR-RL06M v02),記為CSR-M[9]。除橢球效應外,CSR-M采用與球諧系數模型解相同的改正和替換。

1.2 利用球諧系數產品反演全球質量場變化

以等效水柱高(EWH)表示的全球質量場變化Δh可由GRACE/GFO 球諧系數模型計算[2]:

[ΔCnmcos(mλ)+ΔSnmsin(mλ)]

(1)

式中,θ、λ分別表示地心余緯與經度;nmax為截斷階數,在本文中為60;R為地球平均半徑;ρe和ρw分別表示地球平均密度和水密度;ΔCnm和ΔSnm為歸一化球諧系數變化量;Pnm為n階m次歸一化締合勒讓德多項式;kn為n階勒夫數;Wnm為n階m次平滑核函數。

1.3 PnMl方法及其改進方法

GRACE/GFO球諧系數模型在空間域具有明顯的南北條帶誤差[10]。研究表明[3],條帶誤差與球諧系數的次相關誤差模式有關。在以往研究中,為抑制該誤差,Chen等[4]將PnMl(P4M6)方法應用于60階次截斷的GRACE球諧系數模型。

PnMl方法共有3個自由度,Nmin為球諧系數參與擬合部分的最低階數,Nmax為球諧系數參與擬合部分的最高階數,P表示用于擬合相關誤差的多項式次數。數學模型如下:

(2)

(3)

高次位系數的不足使全局保持擬合次數為P的多項式在接近nmax時難以擬合,導致應用60/96階次截斷(nmax=60/96)的GRACE/GFO時變重力場模型在經PnMl方法處理后,大于50/86階次的位系數被完整保留(Nmax=50/86)[11]。未作處理的高階次位系數表現為殘余條帶,嚴重污染經PnMl方法處理的反演結果。改進方法仍保持PnMl方法的3個自由度及其基本思想,對PnMl方法未作處理的高次球諧系數(Nmax,Nmax+1,…,nmax)采取分層擬合。以P4M6方法為例,未經改進時多項式擬合次數P在全局保持不變(P=4),經改進后多項式擬合次數P在不同層次的計算公式為:

(4)

式中,Nmin為PnMl方法中位系數參與擬合部分的最低階數,Nmax為PnMl方法中位系數參與擬合部分的最高階數,nmax為截斷階數。p表示未經改進時在全局保持不變的多項式擬合次數,在改進方法中p=4。

1.4 信噪比評估

采用海陸均方根之比定義信噪比指標[12],計算公式為:

(5)

1.5 基于廣義三角帽的不確定性評估

廣義三角帽方法假設不同觀測序列包含相同的真實信號和相互獨立的噪聲。對任意給定格網點或區域,設Xi(t)(i=1,2,3)為三大機構GRACE/GFO時變重力場模型計算的陸地水儲量變化時間序列,如

(6)

式中,Y(t)為真實地球質量變化信號,σi(t)為不同觀測序列的獨立噪聲。進一步假設任意兩個獨立隨機的觀測值序列之和的方差等于兩者方差之和,構造方差方程:

(7)

若觀測值序列中可用觀測值個數為n,展開式(7)可得n(n-1)/2個方差方程。由于方差方程個數大于或等于觀測數,每個觀測值的噪聲方差可以直接求解,或在n>3時通過最大平方求解。

2 結果與分析

2.1 空域分析

以常用的P4M6方法為例,隨機選取2007-11和2019-09的CSR模型進行測試,得到0.25°×0.25°全球質量變化(以EWH表示)(圖1),圖中Gauss200表示半徑為200 km的高斯濾波。原始時變重力場中存在顯著的南北條帶誤差,經P4M6方法處理后反演結果有所改善(圖1(a)、(d)、(g)),但低緯度區域(30°S～30°N)仍存在殘余條帶。而改進方法可進一步有效抑制殘余條帶誤差,即P4M6方法無法處理的高頻分量中的噪聲(圖1(b)、(h))。由差異圖可知,改進方法可集中處理低緯度區域的條帶誤差(圖1(c)、(i))。對比傳統方法可知,改進方法結合半徑為200 km的高斯濾波可使反演結果更加清晰(圖1(e)、(f))。殘余條帶的減少使后續處理只需應用更弱的濾波策略便可得到清晰的反演結果,使潛在的地球物理信號得以被保留和揭示。

2.2 頻域分析

為分析改進方法在頻域的性能,采用2007-11的CSR模型計算60階次球諧系數階方差(圖2(a)),并給出不同處理方案的振幅及其差異(圖2(b)～(e))。原始高階球諧系數的振幅出現較大異常值(圖2(b)),導致從25階開始,原始階方差呈現異常抬升變化(圖2(a))。經P4M6方法處理后階方差的異常變化被抑制(圖2(c)),但較大的高階信號階方差會導致殘余條帶誤差。而改進方法在50階后的階方差較P4M6方法有所下降,可去除高階突變值(圖2(d)、(e))。

圖2 60階次球諧系數階方差和在不同處理策略下的振幅

2.3 時域分析

在空域和頻域中,改進方法已被證明優于P4M6方法。為進一步比較2種方法在時域的性能,在全球范圍內選取4個經典區域(3°×3°矩形質量塊),矩形中心的大地坐標分別為(44.5°E,21.5°S)、(32.5°E,0.5°S)、(85.5°W,12.5°N)、(0°E,7.5°N),分別位于馬達加斯加西南部、非洲維多利亞湖西北部、尼加拉瓜西北部、非洲沃特湖北部。

按式(1)計算得到矩形區域平均質量變化時間序列,并加入CSR-M作為可靠數據源計算矩形區域質量變化時間序列(圖3)。原始時間序列中的不合理峰值在采用P4M6方法后被部分移除,但仍存在殘余突變值。而改進方法可進一步減少時間序列中的孤立異常值,使信號更符合季節性變化,且與CSR-M計算結果具有較好的一致性。

圖3 3°×3°矩形區域質量變化時間序列

2.4 精度評定

改進方法在空域、頻域和時域的性能提升已被證明。為了量化改進方法去條帶誤差的效果,采用廣義TCH計算三大機構模型以標準差(STD)表示的1°×1°全球不確定性格網評估改進方法,統計3組處理結果在全球及4個矩形區域的平均STD并計算SNR(圖4)。

圖4 不同機構模型精度統計

3種模型在采用改進方法后,全球及4個矩形區域的平均STD較P4M6方法均有所改善。在全球范圍內,CSR、GFZ、JPL三種模型的平均STD分別下降20.50 mm、36.40 mm、19.61 mm(圖4(a)),表明改進方法可以進一步降低時變重力場模型的不確定性。此外,三者的SNR在采用改進方法后較P4M6方法分別從1.62、1.19、1.25增加到1.77、1.35、1.39(圖4(b))。

3 結語

為處理GRACE/GFO時變重力場模型的條帶噪聲,本文提出一種改進的PnMl方法。該方法通過處理PnMl方法忽略的高次球諧系數,進一步抑制條帶噪聲。對比分析改進方法與常用的P4M6方法在空域、頻域和時域的反演結果,同時采用廣義TCH和SNR指標對該方法進行精度評定。主要結論如下:

1)在空域中,改進方法可進一步抑制30°S～30°N區域的殘余條帶誤差,使反演結果更加清晰。

2)在頻域中,改進方法可進一步抑制50階后信號階方差的異常抬升,減少高階系數的階方差突變值。

3)在時域中,改進方法可進一步減少時間序列突變值,使信號更符合季節性變化。

4)經統計,改進方法的SNR、全球和局部區域的平均STD均有明顯改善。其中,CSR、GFZ、JPL三種模型結果在全球范圍的平均STD分別下降20.50 mm、36.40 mm和19.61 mm,SNR從1.62、1.19、1.25增加至1.77、1.35、1.39。這表明,改進方法可通過進一步抑制條帶誤差來降低反演結果的不確定性。