地球重力场是由地球系统质量决定的，反映了地球物质的空间分布[1]．地球内部物质运动、水循环、大气运动等许多自然物理现象会引起地球质量的迁移，使地球重力场处于不断的变动中．地球约71%被海洋所覆盖，研究海洋时变重力对于研究全球气候变化、地幔物质物理化学变化和地壳及岩石圈的变动等具有重要意义．GRACE(Gravity Recovery and Climate Experiment)卫星实现了对地球重力场的重复测量[2]．但GRACE卫星对地球重力场高频信号敏感性较弱，无法恢复较高阶的地球时变重力场模型[3]．GRACE时变重力数据的空间分辨率约为400 km[4]，数据在空间域存在着南北条带误差，须要通过滤波方法予以削弱[5]，会进一步降低GRACE时变重力数据的空间分辨率．随着卫星测高技术不断发展，特别是测高精度和数据分辨率的提高，为海洋重力异常[6]和平均海平面[7]等的研究提供了有效的数据支撑．卫星测高数据中含有丰富的高频信息，国内外学者使用卫星测高技术获得的海面高数据和垂线偏差数据，再利用最小二乘配置(least squares collocation，LSC)、逆斯托克斯公式、拉普拉斯方程和逆Vening-Meinesz(IVM)公式等方法解算静态的海洋重力异常[8-9]．卫星测高数据不仅被用来研究静态的海洋重力异常、平均海平面和海底地形，还被用来研究动态的海平面变化信息[10]，进而研究全球气候变化，但目前还没有利用卫星测高数据反演时变重力的相关研究．CryoSat-2卫星测高数据以其较高的分辨率和测量精度，在建立全球或区域重力场模型[11]中发挥着重大的作用．目前，CryoSat-2卫星已持续观测并收集了十多年的测高数据，为研究海洋时变重力提供了数据条件．孟加拉湾位于亚欧板块和印度洋板块交界处，地质活动比较活跃，研究孟加拉湾及其周边海域的时变重力，有助于该区域板块构造运动和深海扇系统等研究工作的开展．为利用卫星测高数据研究高空间分辨率的海洋时变重力，选用2010年2月—2020年4月间的CryoSat-2卫星测高数据，反演孟加拉湾及其周边海域的时变重力，利用GRACE时变重力数据检验其反演效果，并在周年信号和线性趋势上，对研究海域的CryoSat-2时变重力进行了分析．1 研究区与数据1.1　孟加拉湾及其周边海域选取孟加拉湾及其周边海域作为研究区，地理位置在80°~100°E，0°~23°N之间．研究海域内有安达曼-尼科巴群岛、斯里兰卡岛和苏门答腊岛，海底地形也比较复杂，有印尼海沟、无底峡谷、恒河峡谷、孟加拉盆地、克里希纳-戈达瓦斯盆地和科罗曼德尔盆地等．研究区东部海域位于亚欧板块和印度洋板块交界处[12]，地质构造活动比较活跃．1.2　Cryosat-2卫星测高数据使用AVISO(Archiving，Validation　and Interpretation of Satellite Oceanographic)发布的3.0版本的Cryosat-2沿轨海面高(sea surface height，SSH)数据，该数据存储在采样频率为1 Hz的L2p(Level 2 plus)产品中，数据采用的参考椭球为WGS-84参考椭球．L2p产品数据已针对仪器误差、环境扰动、海况偏差、潮汐影响和大气压力进行了校正[13]．利用CryoSat-2卫星测高数据研究孟加拉湾及其周边海域的时变重力，故将2011年2月到2020年4月间的沿轨SSH数据按月划分成了111组．每个月约有25条上升轨迹和25条下降轨迹通过研究海域，共计约1.4×104个数据．由于执行大地测量任务的CryoSat-2卫星轨道是漂移的，因此每月的观测轨迹并不是有规律的重复，使得不同月份间的测高数据在空间分布上存在着差异．1.3　稳态海面地形数据全球稳态海面地形是卫星测高的关键参考面，测高获得的海面高扣除稳态海面地形值可以得到大地水准面高．CNES(Centre National d'Etudes Spatiales)发布了格网分辨率为7.5′的全球稳态海面地形模型CNES_CLS18．该模型由GOCO05S大地水准面高模型和CNES_CLS15平均海面高模型解算得到[14]，其中GOCO05S大地水准面高模型使用了完整的GOCE任务数据和10.5 a的GRACE任务数据．CNES_CLS18还综合使用了25 a测高数据、温盐数据、浮标数据和水文模型等数据．1.4　参考重力场模型在反演月重力异常时使用了移去-恢复法，因此需要选择一个高精度的参考重力场模型．XGM2019e是一个组合的全球重力场模型，由NGA(National Geospatial-Intelligence Agency)发布，该模型的数据源主要包括GOCO06S重力场模型和NGA编制的地面重力格网数据．一般认为，XGM2019e模型在海洋上的性能更佳，并且独立于现有的高分辨率重力场模型[15]．采用XGM2019e_2159作为参考重力场模型，该模型中球谐系数的阶次完全展开至2 159，阶数扩展至2 190，是目前最新的重力场模型之一．1.5　船载重力异常数据船载测量受海面地形影响较小，通常用于测量区域的重力异常，也被用于检验基于卫星测高数据反演的重力异常．使用NCEI(National Centers for Environment Information)提供的船载重力异常来研究反演的月重力异常是否满足精度要求．孟加拉湾及其周边海域内共有27条船载航线，时间跨度为1963~2008年．由于各船载航线的测量时间、测量仪器和参考椭球不同，船载重力异常数据中存在长波系统误差．选取XGM2019e_2159作为参考重力场，应用二次多项式求解船载数据的改正项[16]，即Δgship=x0+x1Δt+x2Δt2, (1)式中：gship为船载重力异常；Δgship为船载重力异常与XGM2019e_2159的差值；x0，x1和x2为待拟合参数；Δt为测量船观测时间与出发时间的时间间隔．统计船载重力异常改正前后与XGM2019e_2159模型重力异常差值的均方根，结果见表1．10.13245/j.hust.230304.T001表1船载重力异常与XGM2019e_2159模型重力异常差值的均方根测量数据数据集AB测线数量2724测点数量/1045.825 65.3665改正前差值的均方根(mGal)15.1466.481改正后差值的均方根(mGal)14.0825.080将全部船载航线数据记为数据集A，将剔除表2中三条航线的剩余船载航线数据记为数据集B．表1显示：数据集A中的船载重力异常在二次多项式校正后，与XGM2019e_2159模型重力异常差值的均方根比校正前有所减小，但均方根仍较大．数据集B中的船载重力异常进行二次多项式校正后，与XGM2019e_2159模型重力异常差异的均方根降至5.080 mGal，说明表2三条航线中存在着质量较差的观测数据，因此选用数据集B中的船载重力异10.13245/j.hust.230304.T002表2剔除的船载航线的信息v3503so198-1gso198-2g观测时间1978.04.09—1978.04.142008.05.07—2008.06.082008.07.07—008.07.27均方根/mGal改正前差值46.23849.38657.920改正后差值39.74148.46055.117重力仪GSS2-12L&R Seasys S40L&R Seasys S40常数据来检验每个月的重力异常的精度．1.6　GRACE时变重力场美国德克萨斯大学空间研究中心(UTCSR)发布了GRACE和GRACE Follow-On(GRACE-FO) Level-2 RL06 GSM数据产品，基于每月间隔的卫星重力测量数据开发的全球位势模型来估计地球月平均重力场．本研究在国际地球重力场模型中心(International Centre for Global Earth Models，ICGEM)下载了球谐系数阶次展开至60的UTCSR的GSM_DDK3产品，GSM_DDK3数据已扣除了固体潮汐、海洋潮汐、大气潮汐、大气和海洋中的非潮汐变化以及极潮等的影响[17]，并进行了DDK3去相关平滑滤波以削弱条带误差的影响．使用卫星激光测距解算的一阶项、C20与C30对GSM_DDK3数据中相应项进行了替换．利用研究时段内的GSM_DDK3数据计算一个60阶次的平均球谐系数，将GSM_DDK3数据的球谐系数和相应阶次的平均球谐系数相减，得到球谐系数的变化量ΔClm和ΔSlm．根据利用球谐系数确定重力的基本理论[18]得到时变重力函数为      ΔgGRACE(r,λ,θ)=GMr2∑l=0lmaxarl(l-1)×∑m=0lPlm(cosθ)(ΔClmcosmλ+ΔSlmsinmλ), (2)式中：r，λ和θ分别表示计算点的地心向径、经度和余纬；a为WGS-84参考椭球的长半轴；G为牛顿万有引力常数；M为地球质量；l为球谐系数的阶数；m为球谐系数的次数；Plm为完全规格化缔合勒让德函数．利用时变重力函数解算研究海域格网大小为1°的GRACE时变重力数据，该结果中仍包含信号泄露误差和冰川均衡调整(glacial isostatic adjustment，GIA)的影响，因此本研究利用正演模型法(forward-modeling)恢复泄露的信号[19]，利用ICE-6G-C模型[20]对GRACE时变重力数据进行GIA改正．本研究没有对GRACE和GRACE-FO任务之间缺失的11个月数据进行填补．2 研究方法垂线偏差(deflection of vertical，DOV)是解算点重力方向与椭球面法线之间的夹角．计算测高DOV能有效削弱卫星径向轨道误差、海面地形影响以及校正项模型偏差等系统误差，并且DOV中含有丰富的重力场高频成分，有利于恢复高精度高分辨率的海洋重力场，因此可利用DOV数据来反演海洋时变重力．首先，基于沿轨剩余大地水准面梯度利用LSC法解算格网剩余DOV分量；然后，基于格网化剩余DOV分量利用IVM公式计算剩余重力异常．将剩余重力异常与参考重力场模型XGM2019e_2159计算的重力异常相加，得到最终的月重力异常．利用上述方法，反演研究时段内的全部月重力异常数据．对基于CryoSat-2卫星测高数据反演的月重力异常格网数据求均值格网数据并扣除该平均值，进而获得海洋时变重力数据．2.1　数据预处理通过逐月的反演重力异常来研究时变重力．为提高月重力异常的反演精度，可先根据L2p产品手册中定义的高度计、辐射计和地球物理参数的阈值，剔除一些低质量观测数据．同时为降低海平面变化和数据噪声的影响，须对沿轨SSH数据进行高斯低通滤波处理，采用的高斯低通滤波响应函数为h(d)=(2πσ)-1exp[-d2/(2σ2)], (3)式中：d为相邻两点的球面距离；σ为高斯滤波器的宽度．基于全球稳态海面地形模型CNES_CLS18的格网数据，利用双线性内插法推测出观测点上的海面地形值，将滤波处理后的沿轨SSH数据扣除掉海面地形值，得到沿轨大地水准面高N．轨迹上相邻两点P和Q的大地水准面高作差，可有效削弱系统误差，因此计算P和Q两点间的沿轨大地水准面梯度为ePQ=(NQ-NP)/dPQ．(4)将利用卫星测高数据算得的沿轨大地水准面梯度减去XGM2019e_2159模型大地水准面梯度，并剔除两者差异过大的数据，获得沿轨剩余大地水准面梯度，此步骤能够移去重力异常中的长波分量．2.2　利用LSC法计算网格剩余DOV分量LSC是通过计算已知量和估计量的互协方差矩阵以及方差矩阵求估计量的方法．格网点上的剩余DOV的子午圈分量ξres和卯酉圈分量ηres为估计量，格网点位于LSC计算窗口的中心，LSC计算窗口中的沿轨剩余大地水准面梯度eres为已知量．P和Q两点间的剩余DOV分量可以表示为扰动位的函数，因此方差Cξξ(PQ)和Cηη(PQ)可由扰动位的方差解算[21]，进而计算P和Q两点间的协方差Cξe(PQ)和Cηe(PQ)以及方差Cee(PQ)，计算公式为Cξe(PQ)=Cξξ(PQ)cosαPcosαQ;Cηe(PQ)=Cηη(PQ)sinαPsinαQ;    Cee(PQ)=Cξξ(PQ)cosαPcosαQ+ Cηη(PQ)sinαPsinαQ, (5)式中αP和αQ分别为P和Q点处的沿轨方位角．进一步算得LSC计算窗口中各量的互协方差矩阵Cξe和Cηe以及方差矩阵Cee，最后利用LSC计算格网点上的剩余DOV分量：(ξres,ηres)T=(Cξe,Cηe)T(Cee+Cn)-1eres, (6)式中Cn为沿轨剩余大地水准面梯度噪声方差矩阵．2.3　利用IVM公式计算剩余重力异常IVM公式广泛的用于确定海洋重力场[22]．基于格网化剩余DOV分量利用IVM公式计算剩余重力异常，即gres(P)=[GM/(4πR2)]∙    ∬σH'(ψ)(ξres(Q)cosαQP+ηres(Q)sinαQP)dσQ, (7)式中：gres(P)为P点的剩余重力异常；R为地球平均半径；αQP为Q点到P点的方位角；H'(ψ)为核函数的导数，具体为H'(ψ)=-cos(ψ/2)2sin2(ψ/2)+cos(ψ/2)[3+2sin(ψ/2)]2sin(ψ/2)[1+sin(ψ/2)]．(8)因为P和Q两点之间的球面距离ψ无法取0，所以利用IVM计算剩余重力异常时需要考虑内圈带效应的影响，即g0=(s0γ0/2)(ξx+ηy), (9)式中：ξx和ηy为格网剩余DOV子午圈分量、卯酉圈分量的变化速率；s0=(ΔxΔy/π)1/2为内圈带大小，其中Δx和Δy分别为剩余DOV东方向、北方向的格网间隔．将gres和g0相加得到剩余重力异常．3 结果与分析3.1　研究海域时变重力的反演基于Cryosat-2卫星测高数据反演孟加拉湾及其周边海域的时变重力，对沿轨SSH数据进行高斯低通滤波处理，可以有效降低数据噪声的影响．滤波器宽度越大，高斯滤波的平滑程度越好，滤波器宽度过大，可能会剔除非噪声信号，因此必须选择合适的高斯滤波窗口．在没有数据缺失的情况下，沿轨相邻两SSH数据之间的球面距离约为5~7 km，选取0，7和14 km的滤波器宽度进行高斯滤波，统计反演的月重力异常与NCEI船载重力异常差值的最大值、最小值、均值和均方根的平均值，结果见表3．表3显示：当高斯滤波器宽度选取为7 km时，反演的月重力异常与NCEI船载重力异常差值的均方根平均值最小，为5.063 mGal．因此将高斯滤波器宽度固定为7 km，进一步确定其他参数的最优选择．10.13245/j.hust.230304.T003表3不同滤波器宽度反演的月重力异常与船载重力异常差值的统计结果滤波器宽度/km最大值/mGal最小值/mGal均值/mGal均方根/mGal030.153-29.9690.1325.065730.140-29.9650.1315.0631430.071-29.9300.1315.065当计算沿轨剩余大地水准面梯度时，一个月的卫星测高沿轨大地水准面梯度数据较少，须要保证数据的利用率，同时为了提高反演的月重力异常精度，须要剔除与XGM2019e_2159模型大地水准面梯度差异过大的数据．选取不同的沿轨剩余大地水准面梯度剔除阈值，统计反演的月重力异常与NCEI船载重力异常差值的最大值、最小值、均值和均方根的平均值，结果见表4．10.13245/j.hust.230304.T004表4不同剔除阈值反演的月重力异常与船载重力异常差值的统计结果剔除阈值/(″)最大值/mGal最小值/mGal均值/mGal均方根/mGal剔除率/%230.106-29.9340.1285.0678.5330.140-29.9600.1315.0653.5430.140-29.9650.1315.0632.0530.356-29.9640.1325.0701.0由表4可知：当沿轨剩余大地水准面梯度剔除阈值设置为4″时，反演的月重力异常与NCEI船载重力异常差值的均方根的平均值最小，且数据剔除率小于2.0%．当剔除阈值小于3″时，数据剔除率大大增加，数据量减小导致反演的月重力异常精度降低．而当剔除阈值设置为5″时，数据剔除率小于1.0%，沿轨剩余大地水准面梯度中仍存在低质量数据，导致反演的月重力异常精度较低．因此本研究将沿轨剩余大地水准面梯度的剔除阈值设置为4″．利用LSC方法计算格网剩余DOV分量时，选取的LSC计算窗口大小对月重力异常反演结果有很大影响．一个月约50条卫星观测轨迹通过研究海域，相邻两条上升轨迹或相邻两条下降轨迹的间距约为90 km，因此当LSC计算窗口小于50′时，窗口中可能会没有观测数据．选取40′，50′，60′，70′和80′的LSC计算窗口，统计反演的月重力异常与NCEI船载重力异常差值的最大值、最小值、均值和均方根的平均值，结果见表5．10.13245/j.hust.230304.T005表5不同LSC计算窗口反演的月重力异常与船载重力异常差值的统计结果计算窗口大小/(′)最大值/mGal最小值/mGal均值/mGal均方根/mGal4030.065-29.9690.1315.0655030.088-29.9690.1325.0646030.140-29.9650.1315.0637030.154-29.9880.1325.0658030.157-30.0130.1315.066由表5可知：当LSC计算窗口为60′时，反演的月重力异常结果精度最好．当窗口小于60′时，窗口内数据量小甚至没有观测数据，月重力异常结果的精度会降低．当窗口大于60′时，窗口内数据量增多，但窗口内的已知点与待求格网点之间的距离增大，相关性变小，反演的重力异常结果精度也降低．因此本研究选取60′的LSC计算窗口，进一步反演孟加拉湾及其周边海域的时变重力．基于剩余DOV的子午圈分量和卯酉圈分量利用IVM公式计算剩余重力异常时，剩余DOV分量的格网大小会影响剩余重力异常的结果．理论上，IVM公式作为一个积分公式，剩余DOV分量的格网间隔越小，剩余重力异常的计算精度会越高．本研究选取不同格网间隔的剩余DOV分量，统计反演的月重力异常与NCEI船载重力异常差值的最大值、最小值、均值和均方根的平均值，结果见表6．10.13245/j.hust.230304.T006表6不同剩余DOV格网间隔反演的月重力异常与船载重力异常差值的统计结果格网间隔/(′)最大值/mGal最小值/mGal均值/mGal均方根/mGal130.258-30.0320.1335.092230.194-30.0110.1335.082330.140-29.9650.1315.063430.312-29.9100.1375.064530.202-29.9040.1255.104635.245-32.3760.2055.245由表6可知：当剩余 DOV分量的格网间隔为3′时，反演的月重力异常结果精度最好，与NCEI船载重力异常差值的均方根的平均为5.063 mGal．因此基于格网间隔为3′的剩余DOV分量来计算剩余重力异常．最后，恢复XGM2019e_2159模型重力异常，得到格网间隔为3′的孟加拉湾及其周边海域的月重力异常数据．基于划分的111组CryoSat-2卫星测高数据，反演得到111个月重力异常格网数据．利用111个月重力异常格网数据计算一个均值重力异常格网数据，然后将每个月重力异常数据减去这个均值数据，得到研究海域格网间隔为3′的CryoSat-2时变重力数据．由于受GIA影响，海地地壳存在长期形变，因此利用ICE-6G-C模型对CryoSt-2时变重力数据进行了GIA改正．3.2　研究海域时变重力的分析对CryoSat-2和GRACE时变重力数据进行比较，并对研究海域的时变重力进行分析．以2011年2月为例，CryoSat-2月重力场模型和GRACE月重力场模型结果如图1所示，可以发现：CryoSat-2月重力场模型的空间分辨率较高．考虑到近海的卫星测高数据质量较差，因此对距离海岸线小于100 km的近海海域的时变重力数据进行剔除．GSM产品数据进行了DDK3去相关平滑滤波，其光滑效果相当于滤波半径为200 km的高斯滤波处理效果，所以对CryoSat-2时变重力数据进行了半径为200 km的空间高斯平滑处理，以使CryoSat-2时变重力10.13245/j.hust.230304.F001图12011年2月的CryoSat-2月重力场模型和GRACE月重力场模型(色标单位：mGal)数据与GRACE时变重力数据在空间上能够匹配．利用CryoSat-2时变重力数据对孟加拉湾及其周边海域的时变重力进行分析，并和GRACE时变重力数据进行比较，如图2所示．由图2可知：在2011年2月—2020年4月，CryoSat-2时变重力与GRACE时变重力在周年变化上比较相似．利用最小二乘模型来估计研究海域时变重力的周年信号、半年信号和线性趋势。由于半年信号相对较弱，因此这里主要讨论时变重力的周年信号和线性变化趋势，相关统计信息见表7．由表7可知：CryoSat-2时变重力的周年振幅为0.10±0.03 μGal，GRACE时变重力的周年振幅为0.66±0.05 μGal，两者周年振幅具有相同的数量级，周年相位仅偏差约5°，说明两者随时间变化的步调基本一致；研究时段的起点为2011年2月份，CryoSat-2时变重力的周年相位为98.84°±0.43°，说明研究海域时变重力的周年信号约在每年2月份达到最大振幅；CryoSat-2时变重力的线性趋势为0.02±0.01 μGal/a，GRACE时变重力的线性趋势为0.09±0.01 μGal/a，两者基本一致，说明孟加拉湾海域的重力有缓慢增加的趋势．10.13245/j.hust.230304.F002图2孟加拉湾及其周边海域的CryoSat-2时变重力和GRACE时变重力10.13245/j.hust.230304.T007表7研究海域CryoSat-2时变重力信号和GRACE时变重力信号的统计信息时变重力周年振幅/μGal周年相位/(°)线性趋势/(μGal/a)CryoSat-20.10±0.0398.84±0.430.02±0.01GRACE0.66±0.0593.52±0.080.09±0.014 结语基于CryoSat-2沿轨海面高数据解算沿轨大地水准面梯度，将XGM2019e_2159参考重力场模型用于移去-恢复法中，求解格网剩余DOV分量，然后利用IVM公式反演111个月重力异常，最后计算孟加拉湾及其周边海域的3′格网大小的时变重力数据．通过NCEI船载重力异常评估月重力异常的精度，利用GRACE数据检验CryoSat-2时变重力的效果，通过分析得到以下结论．基于CryoSat-2卫星测高数据反演的孟加拉湾及其周边海域3′格网大小的海洋时变重力，在周年信号和长期变化趋势上与GRACE时变重力基本一致，能在一定程度上反映出研究海域2011年2月到2020年4月间的重力变化情况．由于一个月的CryoSat-2卫星测高数据量较少，在剔除沿轨剩余大地水准面梯度的粗差时，剔除阈值设置的过小，数据量会大幅减少，剔除阈值设置的过大，数据中仍然存在着粗差，这都会降低反演的月重力异常精度．因此选择最佳的沿轨剩余大地水准面梯度剔除阈值，可以获得更高精度的月重力异常，有利于海洋时变重力的研究．由于一个月的CryoSat-2卫星测高数据空间分布比较稀疏，在利用沿轨剩余大地水准面梯度计算格网剩余DOV分量时，选取的LSC计算窗口过小，有些窗口内会没有数据，选取的窗口过大，窗口内的已知点与待求点之间的相关性较小，导致月重力异常的精度降低．因此要选择最佳的LSC计算窗口大小，以提高月重力异常的反演精度，进而获得高精度的海洋时变重力数据．以上结果表明：基于CryoSat-2卫星测高数据利用IVM公式能够获得更高空间分辨率的海洋时变重力数据．