卫星重力梯度测量技术的发展，特别是欧空局于2009~2013年实施的地球重力场与静态海洋环流探测卫星计划(GOCE)，为确定高精度高分辨率的地球重力场提供了有效技术支撑[1-2]．GOCE任务是第一颗重力梯度测量卫星，它测定的重力梯度张量更能反映重力场的细部结构和短波变化，有利于建立高精度的静态和局部重力场，因此卫星重力梯度测量数据处理及相关应用研究一直是当前物理大地测量学研究的一个热点问题[3-8]．GOCE任务的科学目标是建立高精度高分辨率的静态重力场模型，其对应100 km分辨率的大地水准面精度为1~2 cm，重力异常精度为1 mGal[2-3]．欧空局GOCE任务下属的高级数据处理部门采用空域法、时域法和直接法研制了多代的GOCE重力场模型[3,9-10]，目前最新发布的第6代模型使用了新版预处理后的GOCE完整观测数据[8]．此外，国内外其他机构利用重力梯度不变量法[5,11]、短弧法[12]和谱组合法[4]等反演了GOCE重力场模型．以上方法均以球谐函数作为基函数反演全球重力场，在高精度局部重力场建模及应用方面存在明显的不足．例如，基于球谐函数建模需要全球覆盖的观测数据，而GOCE观测存在约6.7°的极空白，加剧了重力场反演的病态性[3,11]；同时，病态问题的处理是在全球最优准则下确定统一的正则化参数，往往不能顾及局部重力场的区域特征差异，导致了在信号变化复杂区域出现正则化过度或不足的问题[13-14]．除了以上基于球谐函数的全球重力场建模方法外，国内外学者采用配置法[15-16]、积分法[17-18]和径向基函数法[13-14,19]等，对利用GOCE重力梯度数据建立高精度的局部重力场进行了广泛研究．其中，径向基函数由于具有空域局部集中的特性，在局部重力场建模中具有以下优势：仅须覆盖研究区域的观测数据，解算效率较高；正则化过程中可考虑先验信息，反演精度与可靠性更高；可用于多类不同频段重力观测数据的联合反演．例如，文献[20]模拟研究了不同类型的径向基函数用于局部重力场建模的适用性；文献[13]利用球面样条核函数构建GOCE局部重力场，并提出了依据重力场信号特征的分区正则化策略；文献[14]研究了径向基函数与不同正则化参数选取方法之间的关系，指出Shannon核函数的解算参数与重力场信号强相关；文献[6]利用Poisson小波径向基函数，联合GOCE重力梯度和其他重力数据反演了高精度的局部重力场模型，并评估了重力梯度数据的贡献．此外，卫星重力梯度数据有色噪声的滤波处理也是反演高精度局部重力场的一个关键问题．例如，文献[13]构造完全协方差矩阵抑制有色噪声的影响，但其运算量巨大，计算效率较低；文献[6]和[14]采用带通滤波对卫星重力梯度数据进行预处理，但其丢失了测量频带外的重力场信号，解算过程并不严密．鉴于以上分析，本研究提出以卫星重力梯度数据解算的全球重力场作为背景模型，采用球面样条核函数和分区正则化精化局部重力场的方法及数据处理方案，并且将基于频域内全频段定权的级联滤波方法[21]应用于局部重力场建模．以北美地区GOCE局部重力场反演为例，采用高精度的重力场模型和GPS(全球定位系统)水准数据进行精度评价，验证本文方法的有效性．此外，比较了新、旧两版预处理的GOCE重力梯度数据在重力场反演精度和可靠性的实际效果及差异．1 原理与方法1.1　径向基函数原理将地球表面划分为规则格网，则位于x处的地球外部引力位可表示为在格网中心点位置xi的径向基函数求和形式[13-14,19]V=GMR∑i=1IaiB(x,xi)，(1)式中：GM为地心引力常数；R为地球平均半径；ai为径向基函数系数；I为格网点总数；B(x,xi)为径向基函数，其级数展开形式为B(x,xi)=∑l=0L2l+1R/rl+1blPl(cos ψi)，(2)其中，L为最大展开阶数，r为x的地心向径，Pl为l阶勒让德多项式，ψi为x和xi的球面角距，bl为勒让德系数，它决定了径向基函数的类型．径向基函数格网密度与最大展开阶数L有关，可采用分布均匀的Reuter格网[14]进行计算，并保证全球格网点总数I≥(L+1)2．在局部重力场反演中，由于只须研究区域内的格网点和观测数据参与计算，无法恢复长波重力场信号，因此可采用全球重力场模型作为背景场，从式(1)中扣除后可得T=V-V0=GMR∑i=1IaiB(x,xi)，(3)式中：T为引力位残差；V0为背景场计算的引力位．将式(2)代入式(3)，并对比传统的球谐函数展开式，可得径向基函数系数与球谐函数系数之间的转换关系为[13]flm=∑i=1IaiblY¯lm(θi,λi)，(4)式中：flm为球谐函数系数；Y¯lm为完全规格化的面球谐函数；(θi,λi)为xi位置处的地心余纬和经度．将式(3)在局部指北坐标系(LNRF)下求二阶导数，可得到重力梯度分量的表达式为Tuv=∂2T∂u∂v=GMR∑i=1IaiBuv(x,xi)，(5)式中：Tuv为LNRF下的各重力梯度分量，uv={xx,xy,xz,yy,yz,zz}，即采用GOCE观测精度较高的4个重力梯度分量进行解算；Buv(x,xi)为各梯度分量相应的径向基函数，其具体形式参见文献[19]．须要指出的是，重力梯度数据是在梯度仪坐标系(GRF)下观测得到的，因此须利用坐标转换矩阵[21]将式(5)转换至GRF下建立观测方程进行求解．关于径向基函数的类型，Shannon核函数和球面样条核函数是两类常用的径向基函数，其对应的勒让德系数bl分别为bl=0     (ll0),1     (l0≤l≤L); (6)bl=σl2l+1=Δσl2l+1      (ll0),10-5/l2       (l0≤l≤L), (7)式中：l0为背景模型的最大阶数；σl为重力场信号的阶方差平方根；Δσl为背景模型的误差阶方差平方根．图1给出了Shannon核函数和球面样条(Spline)核函数对应的勒让德系数和空间分布图(具体参数设置见3.1节)．可以看出：这两类核函数的勒让德系数差异较大，但空间分布非常接近，考虑后者采用重力场信号阶方差进行构建，更符合重力场的信号特征，因此采用球面样条核函数进行局部重力场建模．10.13245/j.hust.220920.F001图1Shannon和Spline核函数的勒让德系数及空间分布1.2　有色噪声的去相关滤波由于卫星重力梯度测量误差是具有沿轨强相关性的有色噪声，因此观测值的协方差矩阵为满阵．在最小二乘解算中，通常采用去相关技术构造时域滤波器进行白化滤波处理[22]．常用的高斯-马尔可夫模型为L=Ba+e;DL=σ2P-1, (8)式中：L为观测值向量；B为设计矩阵；a为待估向量；e为误差向量；DL为观测值的协方差矩阵；σ2为方差因子；P为观测值的权阵．将P-1进行LU三角分解，即P-1=RTR，R为正则化矩阵，可采用单位阵代替．并令滤波算子F=(R-1)T，将其对式(8)两边同时进行滤波处理，可得L¯=B¯a+e¯;DL¯=FDLFT=σ2I, (9)式中：L¯=FL；B¯=FB；e¯=Fe．于是，式(9)中的新观测值L¯的权阵为单位阵，即完成了去相关处理．对于等间隔连续观测数据来说，该过程可认为是应用一个离散线性移不变滤波器进行滤波，该滤波器的详细设计见2.2节．1.3　病态问题的正则化方法地球重力场信号的向下延拓及观测值相关等因素，导致利用卫星重力梯度数据反演重力场是一个严重的病态问题．为了获得合理稳定的重力场反演结果，可采用Tikhonov正则化、岭估计和截断奇异值分解等方法进行求解[2,14]．其中，Tikhonov正则化是一种常用且有效的病态问题约束求解技术．根据GRF下各个重力梯度分量去相关滤波处理后的观测方程，采用最小二乘原理和Tikhonov正则化方法可以得到相应求解的法方程为∑uv1σuv2BuvTBuv+1σ2Râ=∑uv1σuv2BuvTTuv, (10)式中：Tuv为GRF下重力梯度分量的观测向量；Buv为各梯度分量观测方程的设计矩阵；â为径向基函数系数的估计向量；σuv2和σ2为方差因子，其中1/σ2也称为正则化参数．根据文献[13]提出的分区正则化策略，正则化矩阵R可根据划分的不同区域i分解为多个正则化矩阵Ri，具体为Rij,j=1      (j∈i);0     (j∉i),式中：j∈i表示â中的第j个系数对应的径向基函数的中心位置在区域i内；j∉i则表示其在区域i以外；且有R1+R2+⋯+Rn=R．为每个区域i单独设置正则化参数，则式(10)可变为∑uv1σuv2BuvTBuv+∑i=1n1σi2Riâ=∑uv1σuv2BuvTTuv，式中σ12,σ22,…,σn2为方差因子，其与径向基函数系数估值可通过方差分量估计迭代求解[23]．若根据重力场信号特征进行分区，可在正则化过程中引入先验信息，进而提高反演精度和可靠性．2 数据处理与背景模型求解2.1　使用数据和处理流程利用2009年11月1日—12月31日共2个月的新、旧两版预处理的GOCE Level 2数据进行重力场反演．首先对EGG_NOM_2和SST_PSO_2数据进行格式转换、粗差剔除和插值等处理，得到重力场解算输入数据，包括重力梯度分量、卫星姿态四元素、简化动力学轨道、惯性系(IRF)与地固系(ERF)坐标转换四元素等．其中：简化动力学轨道用于重力梯度的地理定位；卫星姿态、IRF与ERF坐标转换四元素用于LNRF和GRF坐标系转换．由于重力梯度低频系统误差和有色噪声的影响，利用GOCE观测数据反演重力场模型的低阶位系数精度较低．为此，采用7.5 a GRACE(重力场恢复与气候实验卫星)数据解算的180阶次ITG-GRACE2010S模型法方程联合反演背景重力场模型．以北美地区为例，图2给出了GOCE和GRACE数据联合反演背景重力场，以及基于径向基函数法由背景场和卫星重力梯度数据反演高精度局部重力场模型的数据处理流程．为了验证方法的有效性，采用国际发布的GOCO01S和GOCO06S模型，以及高精度GPS水准数据进行模型精度评价．其中，GOCO01S模型求解采用了与本文同期的GRACE数据和旧版GOCE数据，而GOCO06S模型解算采用了15.5 a的GRACE数据和约4 a的新版GOCE数据，以及Swarm和LAGEOS等卫星重力数据[8]；GPS水准数据密集覆盖了北美地区，且被美国大地测量局(NGS)用于构建GEOID18大地水准面模型(https://geodesy.noaa.gov/GEOID/GEOID18)．10.13245/j.hust.220920.F002图2背景重力场和局部重力场反演的数据处理流程2.2　重力梯度数据滤波处理为了处理GOCE重力梯度数据中的低频系统误差和有色噪声，采用文献[21]提出的移动平均(MA)和自回归移动平均(ARMA)组合的级联滤波器(MA+ARMA)进行处理，其中MA滤波器用于处理低频系统误差，而ARMA滤波器则通过频域内全频段定权的方式抑制有色噪声．MA和ARMA滤波器的差分方程分别表示为c(n)=e(n)-1d+1∑i=0de(n+i)；∑i=0pαic(n-i)=∑j=0qβjw(n-j)，式中：e(n)为输入的重力梯度误差序列；c(n)为MA滤波后的误差序列；w(n)为ARMA滤波后输出的白噪声序列；d为MA滤波窗口长度；p,q为ARMA滤波器的阶数；αi,βj为滤波器系数．与全球重力场反演相比，应用级联滤波器进行局部重力场建模的不同之处在于：研究区域内每个短弧段数据须要单独进行滤波处理，并认为不同短弧段之间的数据不相关；须要对短弧段数据进行扩展，以避免滤波预热阶段舍弃研究区域内的观测值．2.3　背景重力场模型求解全球背景重力场模型不仅作为先验信息用于局部重力场反演，还与采用相同数据反演的局部重力场模型进行比较，以验证径向基函数方法的有效性(见图2)．基于空域最小二乘法建立的重力梯度各分量(Txx,Tyy,Tzz,Txz)法方程与ITG-GRACE2010S模型法方程联合反演背景重力场(详见文献[21])，其中采用MA+ARMA级联滤波器处理重力梯度数据系统误差和有色噪声，采用方差分量估计对各类数据进行定权，并利用Kaula正则化方法进行约束，求解了224阶次的重力场模型位系数．将旧版和新版GOCE数据与ITG-GRACE2010S模型法方程联合反演的背景模型分别记为GOGR-old和GOGR-new，图3给出这两个背景模型和GOCO01S模型相比GOCO06S模型的位系数阶误差(RMS)，表1为相应的全球1°×1°大地水准面和重力异常误差统计．10.13245/j.hust.220920.F003图3相比GOCO06S模型的位系数阶误差10.13245/j.hust.220920.T001表1全球1°×1°大地水准面和重力异常误差统计重力场模型阶次大地准面误差/m重力异常误差/mGal平均值STD平均值STDGOCO01S2240.0030.2080.0466.510GOGR-old2240.0020.2130.0186.612GOGR-new2240.0020.2000.0236.266由图3和表1可知：GOGR-old模型精度略低于GOCO01S模型，但两者在约140阶之后的差异较小，这可能与数据定权和正则化参数选取方法不同有关；而GOGR-new模型优于GOGR-old和GOCO01S模型，并且在约90阶以后的位系数精度有了明显提高，这是由于新版预处理的GOCE重力梯度数据有效削弱了低频误差及粗差的影响，数据质量得到了有效提升．3 北美地区局部重力场反演与分析3.1　局部重力场反演以北美地区[115°W~90°W，35°N~50°N]作为研究区域进行局部重力场建模．为了降低边缘效应的影响，根据试验分析后将GOCE数据和格网点的覆盖范围分别扩展5°和10°．为了降低重力场高频信号混叠的影响，球面样条核函数的最大展开阶数取至240．此外，背景模型的最大阶数取为160，如图1(a)所示，这是考虑到背景场在160阶以后的误差较大，不能反映真实的高频重力场信号．图4(a)为采用GOCO06S模型计算的格网点覆盖范围内的径向重力梯度信号分布，基本反映了研究区域的重力场信号分布；图4(b)为格网点的正则化分区，依据径向重力梯度信号的变化特征划分为4个子区域，各个子区域对应图中格网点的颜色分别为蓝色、黄色、白色和红色．10.13245/j.hust.220920.F004图4径向重力梯度信号分布和格网点正则化分区分别以GOGR-old和GOGR-new为背景模型，采用球面样条径向基函数法和分区正则化策略，反演了北美地区不同系列的局部重力场模型，这些模型解算中采用的背景模型、GOCE数据版本和具体正则化策略如表2所示．10.13245/j.hust.220920.T002表2各局部重力场模型的解算策略重力场模型背景模型GOCE数据正则化GOGR-RBF01GOGR-old旧版未分区GOGR-RBF02GOGR-old旧版分区GOGR-RBF03GOGR-new新版未分区GOGR-RBF04GOGR-new新版分区3.2　高精度重力场模型检核采用高精度的GOCO06S模型和反演的全球与局部重力场模型分别计算研究区域内的格网大地水准面和重力异常，将其差值作为反演模型的大地水准面误差和重力异常误差进行精度评价．图5和图6分别为各种反演模型在北美地区的大地水准面和重力异常误差分布．可以看出：各模型的误差在西部区域较为明显，这是由于西部区域的重力场信号更为复杂，如图4(a)所示；GOCO01S，GOGR-old和GOGR-new在约(102°W，47°N)附近及其东侧的重力场信号相对平滑的区域也存在比较明显的误差，这与全球重力场反演在该区域的正则化不足有关；局部重力场模型GOGR-RBF01和GOGR-RBF03在该区域的误差较小，且分别优于其采用背景模型GOGR-old和GOGR-new，但在西部区域的误差反而有所增大；通过分区正则化处理后的GOGR-RBF02和GOGR-RBF04模型不仅在东部区域的误差较小，在西部区域的误差相比对应的背景模型GOGR-old和GOGR-new也有一定的削弱．10.13245/j.hust.220920.F005图5各模型相比GOCO06S模型计算的北美地区大地水准面误差分布(色标单位：m)10.13245/j.hust.220920.F006图6各模型相比GOCO06S模型计算的北美地区重力异常误差分布(色标单位：mGal)表3为北美地区各模型相比GOCO06S模型的10.13245/j.hust.220920.T003表3北美地区各重力场模型的大地水准面和重力异常误差统计(224阶)重力场模型大地水准面误差/m重力异常误差/mGal最大值最小值RMS最大值最小值RMSGOCO01S0.686-0.6650.16721.850-21.1495.257GOGR-old0.748-0.8220.19223.941-25.8066.023GOGR-new0.609-0.5840.16819.614-18.5055.334GOGR-RBF010.814-0.9330.18525.866-29.4555.702GOGR-RBF020.618-0.5970.16519.148-19.1925.065GOGR-RBF030.738-0.8640.17223.614-27.8735.394GOGR-RBF040.550-0.5180.15417.744-17.2384.839大地水准面和重力异常误差统计．可以看出：分区正则化处理的局部重力场模型GOGR-RBF02和GOGR-RBF04分别优于其全球背景模型GOGR-old和GOGR-new，大地水准面精度分别提升约2.7和1.4 cm，重力异常精度分别提升约1.0和0.5 mGal；当采用相同数据处理方法时，新版GOCE数据反演的重力场模型精度优于旧版数据；与GOCO01S模型相比，同样采用旧版数据反演的GOGR-RBF02模型精度略优，而采用新版数据反演的GOGR-RBF04模型精度提升更为明显．3.3　GPS水准检核为了进一步验证和比较各种反演模型的精度，采用高精度的GPS水准数据对各模型进行外部检核．图7给出了用于检核的3 888个GPS水准点的位置分布，图中显示GPS水准点的分布在东部区域更为密集，但总体上覆盖了大部分的研究区域．表4给出了各种反演模型的GPS水准检核结果，其中平均值表示各模型与GPS水准数据都有一定的系统偏差．10.13245/j.hust.220920.F007图7北美地区GPS水准点位分布(背景为地形)(色标单位：km)10.13245/j.hust.220920.T004表4GPS水准检核结果统计重力场模型阶次平均值STDGOCO06S2240.3640.543GOCO01S2240.3660.566GOGR-old2240.3630.582GOGR-new2240.3690.568GOGR-RBF012240.3740.579GOGR-RBF022240.3690.570GOGR-RBF032240.3720.574GOGR-RBF042240.3690.566m从表4可知：GOCO06S模型的精度为最优，这是由于它采用了更长时间和更多的卫星重力数据进行求解；采用相同的数据处理方法时，新版GOCE数据反演的模型精度优于旧版；分区正则化构建的局部重力场模型GOGR-RBF02和GOGR-RBF04分别优于其背景模型GOGR-old和GOGR-new，且采用旧版数据反演的GOGR-RBF02模型精度提升最为明显；与GOCO01S模型相比，GOGR-RBF02模型的精度略低，而GOGR-RBF04与GOCO01S模型的精度相当，这与3.2节的统计结果有差异，可能与GPS水准点位分布并不均匀有关，因为3.2节是在完全均匀的规则格网点上进行统计的．4 结语以全球重力场模型作为背景场，提出了采用球面样条径向基函数和分区正则化策略由卫星重力梯度数据精化局部重力场的方法和数据处理方案．以2个月的GOCE重力梯度数据和ITG-GRACE2010S模型法方程联合求解的重力场模型为背景场，由GOCE重力梯度数据进一步构建了北美地区高精度的局部重力场模型，精度检核结果表明：采用球面样条核函数和分区正则化方法构建的局部重力场模型能够顾及重力场信号的局部特征差异，解算精度明显优于同期数据反演的全球重力场模型，验证了本文方法的有效性．此外，新版预处理的GOCE数据反演精度明显优于旧版数据，因此建议采用新版数据进行地球重力场建模及应用研究．