海平面的变化是衡量全球气候变化的一个重要指标，与人类活动息息相关．对海平面变化的研究也是地球系统变化研究中非常重要的一个课题，相关的研究在持续进行中，关于海平面变化的具体数值还存在较大争议，因此继续加强海平面变化的监测乃至预测是值得重视的[1-2]．海洋质量变化是引起海平面变化的关键因素之一．传统的验潮站观测、卫星高程监测手段等并不能直接得到海洋质量的信息，而在21世纪初出现的高时空分辨率的卫星重力测量弥补了这一空白，为海洋领域的研究提供了一种全新的手段．2002年3月，重力恢复与气候实验(gravity recovery and climate experiment，GRACE)重力卫星正式发射并投入使用，采集了届时直至2017年的全球时变重力场数据[3]；作为GRACE任务的后继，GRACE-FO(GFO)于2018年5月正式发射并投入使用．到目前为止，GRACE和GFO重力卫星已经提供了长达约20 a的全球时变重力场月解产品，这为全球海洋质量变化的研究提供了关键的数据支撑．然而，在将GRACE时变重力场产品应用于具体问题之前，要经过一系列数据处理过程，这一过程称作后处理．比如，GRACE(-FO) Level-2产品以球谐系数的形式给出了地球重力场的分布，目前最新版本Release 06 (RL06)的空间分辨率约为300 km，但GRACE(-FO)解算的重力场不包含一阶项，因此须用其他手段的观测值作为补充[4]；另外，GRACE(-FO)对C20和C30的解算精度较低，因此须要使用卫星雷达测距(SLR)数据将C20和C30项替换独立于GRACE(-FO)的精度更高的观测值[5-6]．此外，由于给出的产品，尤其在高阶系数中，包含相关性噪声[7]，其在空间分布上表现为南北条带误差，因此在使用GRACE Level-2重力数据解释地表质量变化之前，须要对其做一些相应的滤波处理，同时由于滤波所造成的信号泄漏也需要一定的手段去修复．最后，在使用这些数据解释海洋质量变化之前，还须要做一些相应的模型改正，如GIA模型改正[8]、地震模型改正[9]、海洋潮汐模型[10]和GMAM[11]模型等．显然，上述不同的后处理方法对结果造成的影响并不一致；因此，GRACE(-FO)重力场的不确定度除了由观测手段、测量误差带来的信号不确定度之外，还有上述后处理方法引起的不确定度，二者对于所要研究问题最终结果的影响都是不可忽略的．对于前者而言，由于缺乏其他观测数据的验证，对GRACE以及GFO的噪声水平的估计是比较困难的[12]；然而，对于后者来说，由后处理引起的结果差异，则是可以量化的．已有的一些研究分析了某些后处理步骤对结果的影响，如文献[13]分析了不同的滤波方法对估算全球平均海洋质量变化结果的影响，文献[14]则讨论不同信号泄漏恢复手段对各流域结果的影响．这些研究局限于分析某一特定步骤所用方法对结果的影响大小，而没有对整个后处理流程做全面的分析．本研究将以全球海洋质量变化作为案例，使用不同的后处理手段，得到相应的全球海洋质量变化结果，分析了不同后处理手段对结果的影响，进而给出GRACE的后处理不确定度的量化结果．这种量化结果将明确回答GRACE是否在海洋应用上具有可靠性，为海洋学研究提供重要的参考．1 数据与方法1.1　GRACE卫星重力数据反演全球海洋质量变化1.1.1　重力场数据GRACE Level-2产品数据主要由美国德克萨斯大学空间研究中心(Center for Space Research，CSR)、喷气推进实验室(Jet Propulsion Laboratory，JPL)和德国地学中心(Helmholtz-Centre Potsdam-German Research Centre for Geosciences，GFZ)三家机构提供．除了上述三家机构提供的GRACE时变重力场数据外，也使用了国际时变重力场联合服务(International Combination Service for Time-Variable Gravity Fields，COST-G)产品以作比较．在应用GRACE重力场数据于研究区域之前，还须要对GRACE数据做相应的后处理，以及必要的模型修正，图1给出了大致的后处理流程，其中具体的步骤将在后几个小节介绍．10.13245/j.hust.230311.F001图1GRACE Level-2 重力场数据简化处理流程在重力场低阶项数据的替换步骤中，使用官方提供的TN-13文档的一阶项以及TN-14的C20 项和C30项作为GRACE低阶项的替换．1.1.2　滤波方法Swenson等发现了GRACE产品在高阶系数的奇次项之间和偶次项之间存在相关误差，其在空间域则表现为“南北条带误差”，并提出了相应的去相关滤波方法[7]，该方法将特定次不同阶的系数作多项式拟合，并从原始系数中扣除．Kusche等[15-16]提出的DDK滤波方法也可以很好地减弱这样的误差．去相关滤波并不能很好地处理高纬度地区的“条带”，还须进一步引入高斯滤波器[17]、各向异性高斯滤波器[18]、扇形滤波器[19]等低通滤波器来进一步抑制噪声．使用上述不同种类和参数的去相关滤波和低通滤波组合来抑制高阶系数的相关误差并比较其差异．1.1.3　信号泄漏的处理不同滤波手段的引入虽然很好地抑制了高阶噪声，但不可避免地对高阶信号部分也有所削弱，在空间上则表现为信号的偏差(bias)和泄漏(leakage)，前者指研究区域内的信号往外流失，后者则是指外部区域的信号流至内部．不同区域信号的偏差和泄漏影响程度不一．对于偏差的修复，一般可以根据不同手段预估出一个尺度因子来修正信号；对于泄漏而言，文献[17]提出了一种简单的处理手法(这里称二次滤波法)来估计并扣除．此外，文献[20]也提出了正向建模法来处理信号泄漏误差．也有研究在研究陆地水文信号时，并不区分偏差和泄漏，而是将二者统称为信号泄漏来统一处理，模式数据估计区域或网格化尺度因子[21]．特别地，在处理全球尺度的海洋信号时，由于研究区域尺度较大，文献[10]在求和海洋信号时，对海洋区域的选取直接忽略了距海岸线300 km的区域来回避陆海交界处的信号泄漏．对于信号泄漏的恢复策略，也将使用上述几种不同的方法来恢复海洋信号的偏差和泄漏误差，并比较其差异．1.1.4　模型修正冰川均衡调整(Glacial Isostatic Adjustment，GIA)是指近十万年来，由于地球对冰盖与海洋之间的质量再分配的黏弹性响应，地球重力场主要受长期趋势的影响．由于全球海洋尺度较大，其结果受到GIA影响将尤其明显．表1给出了GIA效应对GRACE结果在全球海洋范围内的影响，其中，扣除300 km缓冲带的海洋区域由文献[10]给出．在用GRACE时变重力场数据解释海洋质量变化之前，须扣除这一影响．10.13245/j.hust.230311.T001表1不同GIA模型对GRACE观测的海洋质量趋势影响模型全球海洋区域扣除300 km缓冲带的海洋区域A13[22]ICE-6G_C[23]ICE-6G_D[24]Caron18[8]-0.93-0.84-0.87-1.00-1.16-1.10-1.11-1.36mm/a当解算GRACE二级产品GSM模型时，大气和海洋信号作为AOD1B模型被扣除．当使用GRACE产品计算海洋质量的时空分布时须要将这一部分信号恢复到模型中．这一部分信号可以通过随GSM产品附带的GAD模型得出；同时，为了与后续测高、Argo数据保持一致，在恢复海洋模型的同时还须做出相应的IB改正[10]．重力场GSM在解算过程中已扣除了AOD1B的信号，此时GSM信号在全球范围内并不守恒．而GSM并没有解算零阶项而将其认为是一个常数(即全球质量守恒)，这实际上引入了GMAM(global mean atmosphere mass)误差．GMAM可以通过随GSM产品发布的GAA产品得出[11]，此项误差也须要从GSM场中扣除．地震可以在短时间内引起大范围的质量迁移，研究海洋质量变化时这一信号须要被扣除．现在一般认为，GRACE可以观测到8级以上的地震信息．Tang等[9]在利用GRACE数据分析全球海洋质量变化之前，利用最小二乘的方法对GRACE观测信号拟合出某一点的地震信息并扣除．参照文献[9]的方法，分析了GRACE期间强于里氏8.5级的地震，扣除了相应的拟合信号，具体为yseis(t)=c0          (tt1);ci+pi(1-e-(1-ti)/τi)          (ti≤t≤ti+1tn);cn+pn(1-e-(1-tn)/τn)          (t≥tn),式中：yseist为对所研究的区域拟合出的地震信号；n为该区域受到地震影响的总次数；ti为第i次地震发生的时间；τi为第i次地震的弛豫时间；ci和pi为须要拟合得出的参数．表2给出了所分析的地震事件的发生时间．10.13245/j.hust.230311.T002表22005~2020年期间陆海交界处的大地震地区发生时间里氏震级Nias，印度尼西亚Maule，智利Tohoku-Oki，日本Sumatra，印度尼西亚鄂霍茨克海2005-03-282010-02-272011-03-112012-04-112013-05-248.68.89.08.68.31.2　联合卫星测高、海水温盐数据反演全球海洋质量变化使用哥白尼海洋环境监测系统(Copernicus Marine Environment Monitoring Service，CMEMS)提供的每月全球海平面异常数据(Sea Surface Hight，SSH)得到全球海平面的变化趋势，SSH融合了Jason-3，Sentinel-3A和HY-2A等多颗卫星的测高数据．使用由中国Argo实时资料中心(CARDC)提供的“BOA_Argo”网格资料[25]，并根据海平面方程计算得到全球海洋的比容海平面变化．SSH数据集中，里海被纳入全球海洋区域．使用SSH数据集时，将首先扣除这一区域的信号，并扣除相应的GIA影响，这一影响在海洋区域的数值约为-0.3 mm/a[11]．海平面高度的变化可以细化为海洋质量的变化和比容海平面的变化，卫星测高系统可以得到海平面的高度变化，而“BOA_Argo”网格资料所提供的海水温盐数据则可以得到其中比容海平面变化的部分．因此，将卫星测高得到的海平面变化扣除比容海平面，可以得到海洋质量变化，这一结果理论上应该与GRACE所观测到的海洋质量变化保持一致．2 方法验证2.1　联合卫星测高和海水温盐数据得到的全球海洋质量变化图2给出了海平面高度变化和比容海平面变化和联合卫星测高数据和海水温盐数据得到的海洋质量变化，表3列出了其数值结果．由于CMEMS提供的卫星测高月解不包含2020年5月份以后的数据，当计算测高和Argo数据时，给出的是2005—2019年间的结果．10.13245/j.hust.230311.F002图22005—2019年间，使用测高卫星和海洋温盐数据估计的全球海平面和海洋质量变化趋势10.13245/j.hust.230311.T003表3联合卫星测高和Argo数据得到的海平面变化海平面和海水质量变化线性趋势/(mm/a)周年振幅/mm周年相位/(°)卫星测高海平面Argo比容海平面测高-Argo3.850.892.965.894.129.37144.506.35161.642.2　重力场低阶项的影响图3给出了2005—2020年期间，COST-G给出的重力场在不同低阶项替换的情况下前三阶对全球海洋质量变化的贡献．表4分别给出了不同机构以及GSFC SLR估计的C20和C30项对于海洋质量等效水高的贡献，其中GSFC给出的结果在2012年之前没有C20的值，这一段时间内使用CSR机构的数据代替．10.13245/j.hust.230311.F003图3使用COST-G数据并做不同的低阶项替换后，前三阶对全求海洋质量变化的贡献10.13245/j.hust.230311.T004表4C20和C30对海洋质量变化等效水高趋势的影响 mm/a数据来源C20项贡献C30项贡献/10-2COST-G0.176-2.65CSR0.211-2.77GFZ0.116-2.60JPL0.141-2.49GSFC0.120-2.622.3　滤波方法和信号泄漏恢复策略的影响图4和表5给出了不同滤波和信号泄漏恢复方法在海洋区域上质量变化的结果．表5中滤波1~滤波4分别为300 km各向同性高斯滤波、500 km各向同性高斯滤波、(300，800) km各向异性高斯滤波和(300，300) km扇形滤波，泄漏恢复a~c分别是正向建模法、使用扣除300 km缓冲带的海洋区域和使用文献[17]提出的信号泄漏恢复方法并乘以相应的尺度因子．图4则给出了其中部分组合策略的时间序列．10.13245/j.hust.230311.F004图4选取不同滤波和信号泄漏恢复策略得到的全球海洋质量变化趋势10.13245/j.hust.230311.T005表5滤波与信号泄漏恢复策略对海洋质量等效水高变化的影响策略泄漏恢复a泄漏恢复b泄漏恢复c滤波1滤波2滤波3滤波41.8911.831.841.862.182.0222.152.171.781.791.7931.78注：上标1表示图4中的组合策略1；上标2表示图4中的组合策略2；上标3表示图4中的组合策略3．mm/a2.4　模型修正的影响对于全球海洋而言，从表1可以看出：GIA的影响是尤为重要的一点．然而，各GIA模型的误差来自于GIA模型本身，因此所有案例均使用文献[24]给出的GIA模型．GMAM信号如图5所示，其代表了解算GSM模型时将0阶项置零时引入的误差，这一部分误差可以从GAA产品导出．实际分析时，将GRACE-GSM海洋信号扣除这一部分信号．10.13245/j.hust.230311.F005图52005—2020年间，使用GAA模型零阶项计算的GMAM信号计算海洋质量变化时，陆海交界处的地震信号也是不可忽视的系统误差，这一部分信号可以由GRACE观测信号拟合而得．以2012年苏门答腊地震为例，图6给出了地震范围内一点(2.5°N，93.5°E)的GRACE观测信号以及相应的拟合信号，实际应用中，将表2所示的地震影响范围内的GRACE信号扣除拟合而得地震信号．表6给出了2005—2020年间地震信号对海洋质量变化的趋势影响．10.13245/j.hust.230311.F006图62005—2020年间，苏门答腊区域一点(2.5°N，93.5°E)的重力场信号及拟合出的2012年地震信号10.13245/j.hust.230311.T006表62005—2020年期间地震对海洋质量变化的影响地震全球海洋扣除300 km缓冲带的海洋区域Nias 2005Maule 2010Tohoku-Oki 2011Sumatra 2012鄂霍茨克海 2013地震总效应0.0090.0090.0220.0290.0130.0950.0060.0050.0250.0180.0090.067mm/a3 不确定度分析第2节所述的每一步策略对结果都有不可忽略的影响．由于后处理流程的方法前后相互影响，每一步骤对最终结果的影响并不是线性的．除了取决于自身的参数外，还取决于其他步骤策略以及对应参数的选择．对GRACE数据处理的每一步骤，选取如表7所示的方法和参数，其中，模型修正则是“有”或者“无”．因此，最终形成的组合数量为4×3×4×3 ×3×3=1 296种，以上多组后处理流程将会相互比较得到结果．最终的结果由图7~图10和表8给出，其中为了与测高-Argo结果比较，表8里的内容均是2005—2019年间的相应结果．10.13245/j.hust.230311.T007表7使用的后处理方法和参数数据选择低阶项替换高斯滤波信号泄漏恢复模型修正去相关滤波COST-G仅替换1阶项300 km高斯滤波正向建模GIA：ICE-6G_D (固定)不使用CSR替换1阶项、C20500 km高斯滤波次滤波法GAD模型P3M5GFZ替换1阶项、C20和C30(300，800) km各向异性陆海缓冲区法GMAM模型滑动窗口拟合JPL—(300，300) km扇形—地震模型固定窗口拟合10.13245/j.hust.230311.F007图7GRACE后处理不确定度10.13245/j.hust.230311.F008图8剔除最大、最小的5%线性趋势后得到海洋质量变化趋势10.13245/j.hust.230311.F009图9GRACE后处理不确定度结果10.13245/j.hust.230311.F010图10不同方法计算全球海洋质量变化趋势频次图(90%样本)10.13245/j.hust.230311.T008表82005—2019年间，GRACE后处理最大、最小趋势(90%样本)海平面和海水质量变化线性趋势(mm/a)周年振幅/mm周年相位/(°)测高-ArgoGRACE最大趋势GRACE最小趋势GRACE平均趋势GRACE中位趋势不确定度2.962.301.661.941.930.649.3710.787.919.549.771.48161.64175.82173.04174.53174.622.78图7为这1 296个案例给出的全球海洋质量在研究时间段内的结果，其中：黑色实线是测高-Argo得到的海洋质量变化结果；红、绿实线是所有组合得到全球海洋质量变化趋势的上、下包络线；中间虚线则是触及到包络线的组合，其余组合的结果则不在图中显示；触及到包络线的部分结果的时间序列以虚线给出，其他部分的结果则不在图里给出．同时，为了规避一些少数极端的结果，图8给出了类似于图7，但根据最终计算的趋势大小，扣除了前后各5%的样本的结果，其中GRACE结果仅包含上下包络线．从图7、图10和表8可以看出：GRACE在海洋上的结果与测高-Argo的结果总体基本一致，但是在细节上存在一些不符合的现象，如相对于符合较好的周年振幅而言，GRACE与Argo-测高的结果在线性趋势上并十分相符．一方面从数据本身所代表的物理信息而言，Argo数据本身不包含2 000 m水深以下的信息，同时，Argo浮标也并没有对全球范围内的海洋区域都覆盖到[25]；并且通常在海洋区域，扣除了已有的模型后和进行一些改正之后，才认为得到的信号是由海洋质量变化引起的，但是这些模型和相应的改正本身存在较大的不确定度．图9给出了图8中上包络线与下包络线的差值及其拟合的线性趋势、年周期结果，其中：黑色实线是图8中上下包络线的差值；红色虚线则是其拟合的线性趋势和年周期信号，后面称这个差值为“后处理不确定度”．从图9的结果可以看出：在研究的这段时间内，后处理不确定度随时间有增大的趋势，尤其在GFO期间，后处理不确定度表现得比GRACE期间更为明显，这可能与2016年以后GRACE和GRACE-FO期间的单加速度计模式[26]有关．由表8可见：通过计算90%样本的最大、最小趋势的差，可以得到后处理得到全球海洋质量增长趋势不确定度约为0.64 mm/a．另外，经统计计算，所计算的样本中，线性趋势的标准差约为0.24 mm/a．除了线性趋势之外，图9还显示出后处理不确定度具有年周期性．从拟合结果可见：其振幅为1.48 mm，并且在每年的3月份左右，GRACE的后处理不确定度较小，而9月份左右则达到最大峰值．这表示在每年3月份左右GRACE在海洋上的结果相对更为可信．除了上述对时间序列的分析之外，图10还给出了上述1 296种的结果所得到的线性趋势在扣除最大、最小的5%样本后的频次图．通常而言：期望这样的统计结果呈现高斯分布是因为这样能够将后处理不确定度解释为随机误差．但是，从图10可以看出：对于不同组合得到的线性趋势呈双峰分布，这说明这里面有一些主导的关键物理机制．经统计分析发现：这样的双峰差异主要来自于后处理的信号泄漏恢复步骤中，是否使用陆海缓冲带的方法，因此图10中的结果进一步以使用缓冲带和非缓冲带的信号泄漏恢复方法加以区分．从表1可以看出：使用陆海缓冲带后得到的GIA趋势比整个海洋区域计算的GIA区域更大．因此可以认为，这样的双峰结果很大一方面来自于GIA模型在不同定义的海洋区域的结果差异．实际上，即使考虑了GIA模型在全球海洋区域和300 km缓冲带的海洋区域的差异，这样的双峰现象虽然得以改善，然而仍然是存在的，从这一角度来看，现有的GRACE后处理流程存在系统误差，如何改善或提出更有效可行的后处理方法，是值得关注的一个话题．4 结论以全球海洋质量变化为案例，使用GRACE和GRACE-FO时变重力场数据，通过不同的后处理方法组合得到了GRACE后处理引起的海洋质量增长趋势和年周期振幅的最大、最小值和不确定度．GRACE在海洋上的线性后处理不确定度约为0.64 mm/a，标准差约为0.24 mm/a，并且其周期性显示，GRACE在每年3月的可信度相对优于9月份的可信度．总体来看：GRACE在海洋应用上的后处理不确定度在一个较稳定的范围，因此认为GRACE的重力场数据是可靠的、有意义的．从结果也可以看出：GRACE在海洋上的结果与测高-Argo结果在2015年之前较为相符，后续时间则符合得并不理想，这一方面由于Argo数据和GRACE数据各自在测量原理上不能完美地表示整个海洋的信号或将海洋信号与其他信号(如固体地球信号等)分离开，另一方面这与GIA模型的不确定性也有很大的关系．通过以上研究可以看出：GRACE后处理不确定度是不可忽略的，并且呈现一定的规律，说明这些处理方法存在系统性误差，因此对后处理研究还须要进一步深入，才能得到更加准确的GRACE反演结果．除了海洋的研究之外，GRACE在其他小尺度流域等的研究可能存在更大的后处理误差，这将会对一系列的科学应用带来更大的挑战．GRACE后处理不确定度的研究对评估GRACE反演地表水质量结果有着重要的意义，首次在全球海洋区域做了一定的分析，后续将进行更为广泛细致的研究．