联合GRACE和GNSS的局部地表质量变化反演方法及系统

摘要

本发明提供一种联合GRACE和GNSS的局部地表质量变化反演方法及系统，包括获取由局部地表质量变化引起的GRACE星间重力位差和GNSS垂直位移观测值；根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程，得到相应法方程；根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程，得到相应法方程；根据最小二乘联合平差反演局部地表质量变化，包括根据法方程结合先验信息方程形成GRACE和GNSS联合反演模型，然后给定两类观测值噪声方差以及正则化参数初始值，利用方差分量估计通过迭代计算确定最优权比，得到联合反演的局部地表质量变化。本发明提高了局部地表质量变化反演的精细度和可靠性。

说明书

技术领域

本发明涉及联合空间大地测量观测数据反演局部地表质量变化的方法，特别涉及一种联合GRACE星间重力位差和GNSS垂直形变数据反演精细局部地表质量的方法及系统。

背景技术

地表质量变化主要包括陆地水储量变化、极地冰盖和山岳冰川融化、大气压变化、海洋质量变化以及和固体地球物理现象相关的其他质量变化现象。这些全球或局部的地表质量迁移与重新分布，不仅会引起地球重力场的改变，还会导致固体地球表面产生形变。因此，通过观测地球重力场和地球形变场的时空变化可以推演和监测地球系统的物质迁移和交换过程，并且获取时变重力场与形变场的精细程度越高，其包含的地球系统物质迁移信息则越丰富，这对大地测量学、固体地球物理学、地球动力学、海洋学、冰川学以及全球环境变化等研究都具有十分重要的物理意义

近年来，卫星重力测量技术的实现，特别是美德联合实施的GRACE计划，为准实时获取全球高精度时变重力场进而反演地表质量变化提供了新途径。不同于传统的重力探测技术，GRACE能够获取全球高精度、均匀覆盖的月时变重力场，可在季节时间尺度上以1cm等效水高精度、400km空间分辨率提供地表质量变化信息

由于地球质量重新分布导致的地表形变能够被GNSS以毫米级精度连续地观测，因此可以利用GNSS负荷形变资料独立地反演地表质量变化

目前国内外学者的研究成果侧重于单独使用GRACE或GNSS数据反演局部地表质量变化，或对两类数据的反演结果进行一致性和相关性分析。虽然Adusumilli等(2019)在反演美国周边地区水储量变化时利用GRACE和GNSS观测资料进行联合反演，但其实质是使用GRACE反演结果对GNSS反演结果进行空间约束，并没有从观测数据层面进行联合反演

本发明涉及的参考文献如下：

[1]Blewitt G.,D.Lavallée,P.Clarke,K.Nurutdinov.A new global mode ofearth deformation:seasonal cycle detected.Science,2001,294(5550):2342-2345.

[2]Tapley,B.D.；Bettadpur,S.；Ries,J.C.；Thompson,P.F.；Watkins,M.M.GRACEMeasurements of Mass Variability in the Earth System.Science 2004,305,503–505.

[3]许厚泽,陆洋,钟敏,郑伟,张子占.卫星重力测量及其在地球物理环境变化监测中的应用.中国科学:地球科学,2012,42(6):843-853.

[4]Swenson,S.；Wahr,J.Post-processing removal of correlated errors inGRACE data.Geophys.Res.Lett.2006,33,L08402.

[5]So′snica,K.；Jaggi,A.；Meyer,U.；Thaller,D.；Beutler,G.；Arnold,D.；Dach,R.Time variable Earth’s gravity field from SLR satellites.J.Geod.2015,89,945–960.

[6]Wu,X.Large-scale global surface mass variations inferred from GPSmeasurements of load-induced deformation[J].Geophysical Research Letters,2003,30(14):253-266.

[7]Zhong B,Li X,Chen J,et al.Surface Mass Variations from GPS andGRACE/GFO:A Case Study in Southwest China[J].Remote Sensing,2020,12(11):1835.

[8]Argus,D.F.；Fu,Y.；Landerer,F.Seasonal variation in total waterstorage in California inferred from GPS observations of vertical landmotion.Geophys.Res.Lett.2014,41,1971–1980.

[9]Fu Y,Argus D F,Freymueller J T,et al.Horizontal motion in elasticresponse to seasonal loading of rain water in the Amazon Basin and monsoonwater in Southeast Asia observed by GPS and inferred from GRACE[J].Geophysical Research Letters,2013,40(23):6048-6053.

[10]Adusumilli S,Borsa A A,Fish M A,et al.A Decade of Water StorageChanges Across the Contiguous United States From GPS and Satellite Gravity[J].Geophysical Research Letters,2019,46(22):13006-13015.

[11]Han S C,Shum C K,Braun A.High-resolution continental waterstorage recovery from low–low satellite-to-satellite tracking[J].Journal ofGeodynamics,2005,39(1):11-28.

[12]Zhong B,Li Q,Chen J,et al.Improved Estimation of Regional SurfaceMass Variations from GRACE Intersatellite Geopotential Differences Using aPriori Constraints[J].Remote Sensing,2020,12(16):2553.

发明内容

针对现有单一空间大地测量观测数据反演局部地表质量变化方法的不足，本发明提供了一种联合GRACE星间重力位差和GNSS垂直形变数据反演精细局部地表质量的方案。

本发明采用的技术方案为一种联合GRACE和GNSS的局部地表质量变化反演方法，包括如下步骤，

步骤1，获取由局部地表质量变化引起的GRACE星间重力位差和GNSS垂直位移观测值；

步骤2，根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程，得到相应法方程；根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程，得到相应法方程；

步骤3，利用步骤2所得结果根据最小二乘联合平差反演局部地表质量变化，包括根据法方程结合先验信息方程形成GRACE和GNSS联合反演模型，然后给定两类观测值噪声方差以及正则化参数初始值，利用方差分量估计通过迭代计算确定最优权比，得到联合反演的局部地表质量变化。

而且，步骤2中，

根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程如下，

其中，y

根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程如下，

其中，y

而且，步骤3中，设增加先验约束方程

x

其中，x

采用的联合反演模型为，

其中，P

给定噪声方差

而且，根据步骤3中联合反演模型反演局部地表质量变化，与单独利用GRACE或GNSS反演结果进行比较，以反演结果与真实信号残差的标准差以及阶方差RMS作为反演结果精度的评判。

另一方面，本发明还提供一种联合GRACE和GNSS的局部地表质量变化反演系统，用于实现如上所述的一种联合GRACE和GNSS的局部地表质量变化反演方法。

而且，包括以下模块，

第一模块，用于获取由局部地表质量变化引起的GRACE星间重力位差和GNSS垂直位移观测值；

第二模块，用于根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程，得到相应法方程；根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程，得到相应法方程；

第三模块，用于利用第二模块所得结果根据最小二乘联合平差反演局部地表质量变化，包括根据法方程结合先验信息方程形成GRACE和GNSS联合反演模型，然后给定两类观测值噪声方差以及正则化参数初始值，利用方差分量估计通过迭代计算确定最优权比，得到联合反演的局部地表质量变化。

或者，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种联合GRACE和GNSS的局部地表质量变化反演方法。

或者，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种联合GRACE和GNSS的局部地表质量变化反演方法。

与现有技术相比本发明具有以下优点和有益效果：

1、提高了局部地表质量变化反演的精细度和可靠性。

受限于GNSS测站分布或GRACE频谱敏感范围有限等因素，当前单独利用GRACE或GNSS数据反演局部地表质量变化都存在一定的缺陷与不足。本发明联合GRACE和GNSS反演局部地表质量变化的方法，可以充分地利用两类观测数据的优点，实现优势互补，进而提高了反演结果的精细度和可靠性。

2、方案的理论依据充分，实现稳定有效。

通过分别对GRACE、GNSS反演结果以及GRACE和GNSS联合反演结果与原始信号的残差标准差以及阶方差RMS进行对比可知，联合反演结果与原始信号残差更小，并且在高阶部分阶方差RMS的信噪比有了明显提高，反映出具有更高的空间分辨率。本发明充分利用GRACE和GNSS数据的各自优点，因而联合反演方法对提高局部质量变化的精细度和可靠性是有效的。

附图说明

图1为本发明实施例的方法流程图。

具体实施方式

以下结合附图和实施例对本发明技术方案进行具体描述。

本发明充分考虑GRACE和GNSS技术反演局部地表质量变化的优缺点，实现优势互补，通过联合GRACE星间重力位差和GNSS垂直位移观测数据反演精细的局部地表质量变化。由此本发明提出了一种联合GRACE星间重力位差和GNSS垂直形变数据反演精细局部地表质量的方法。

下面以南美洲作为实验区域，选取2005年9月的GLDAS水文模型计算的水储量变化作为原始质量变化信号，通过闭合数值模拟并结合图1说明本发明实施例提供的一种联合GRACE和GNSS的局部地表质量变化反演方法，具体流程步骤实现如下：

步骤1，获取由局部地表质量变化引起的GRACE星间重力位差和GNSS垂直位移观测值。

根据牛顿万有引力定律，GRACE双星中点处的星间重力位差与局部地表质量变化之间的关系可以表示为：

其中，

T

G为万有引力常数：6.673×10

k

N为研究区格网点数，j为格网点标号，j＝1,2,…N；

根据质量负荷的格林函数理论，GNSS垂直位移与局部地表质量变化之间的关系可以表示为：

其中，u为GNSS垂直位移；M

实施例中，为了方便计算，研究区域被划分为2°×2°的均匀格网，包括海洋区域在内共有1216个格网。以2005年9月的GLDAS水文模型计算的水储量变化作为原始信号，根据研究区的经纬度范围和GRACE卫星星历数据计算得到GRACE星间重力位差数据，包括海洋区域共有38696个观测值。通过方程

步骤2，根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程；根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程；根据GRACE和GNSS观测方程形成联合反演法方程组。

GRACE星间重力位差的观测方程为：

其中，y

根据GRACE星间重力位差的观测方程便可得到其法方程：

其中，

GNSS垂直位移的观测方程为：

其中，y

由于利用GRACE或GNSS数据反演局部地表质量变化属于离散病态问题，可利用先验信息对病态方程进行约束求解。根据GRACE或GNSS的观测方程，增加先验约束方程

x

其中，x

根据GNSS垂直位移的观测方程便可得到其法方程：

其中，

单独使用GRACE或GNSS反演局部地表质量变化的反演方程

||Ax-y||

其中，α为正则化参数，A为设计矩阵，y表示观测值，min表示最小(具体可理解为当目标函数最小时则可以求出待估的参数x)。

待求的地表质量变化

这里P为权阵，通过给定噪声方差

实施例中，GRACE星间重力位差的观测方程为：y

步骤3，通过步骤2所得GRACE和GNSS观测方程的法方程组根据最小二乘联合平差反演局部地表质量变化。

GRACE为“物理”观测量，GNSS为“几何”观测量，为实现这两类观测值的最佳融合，需要充分顾及它们在时空分辨率、空间覆盖和频谱敏感性等方面的差异，合理考虑其最优融合的模式及准则，本发明通过方差分量估计确定两类数据的最优权比。

联合GRACE和GNSS反演局部地表质量变化属于离散病态问题，需要利用正则化方法对病态法方程组进行约束求解，常用的包括Tikhonov正则化、截断奇异值分解(TSVD)以及岭估计等方法都受限于如何快速有效地获取最优正则化参数，本发明以先验地球物理模型信息(如水文模型)作为约束，通过最小二乘迭代方法根据观测数据本身自适应地、有效地获取正则化参数。

本发明进一步提出，联合反演需要确定联合反演模型和定权，还需包含以下子步骤：

3.1根据法方程并结合先验信息方程形成GRACE和GNSS联合反演模型，即形成联合反演目标函数：||A

则待求的地表质量变化参数估计值的联合反演模型可以表示为：

3.2给定两类观测值噪声方差以及正则化参数初始值，利用方差分量估计通过迭代计算确定最优权比，可以计算出联合反演的局部地表质量变化。

给定噪声方差

本实施例中，分别利用GRACE、GNSS以及联合GRACE和GNSS反演结果与原始信号的残差标准差和阶方差RMS及信噪比进行统计分析便可以评估联合反演结果精度。

实施例中，联合GRACE和GNSS含误差模拟数据反演南美洲局部地表质量变化。具体包含以下子步骤：

3.1根据步骤2中GRACE与GNSS的观测测方程及先验约束方程则可以得到联合反演反演模型：

3.2以模拟的GRACE和GNSS噪声标准差作为方差初始值并给定一个正则化参数，通过方差分量估计迭代定权便可以估计最优的局部地表质量变化结果。得到GRACE和GNSS联合反演结果，由于GNSS与GRACE在数据分布以及频谱敏感范围等方面各有优势，相对于GNSS反演结果而言，在没有测站分布区域，联合反演结果中GRACE的贡献很好地补充了该区域的地表质量变化信息；相对于GRACE反演结果而言，在有GNSS测站分布的区域(特别是GNSS密集分布)，GNSS的贡献很好地恢复了局部高频信号，提高了反演结果的精细度。可得到为联合反演结果与原始信号的残差图，其中反演结果与原始信号残差标准差为29.40mm。可得到GNSS、GRACE以及二者联合反演的结果与原始信号残差的阶方差RMS，根据结果可知，GRACE反演结果在低阶项(大约30阶之前)优于GNSS反演结果，而GNSS反演结果大约在30阶以后优于GRACE反演结果，将两者联合反演的结果大约在15阶以后优于GRACE或GNSS独立反演结果。同时，从阶方差RMS的信噪比可以看出，GRACE和GNSS联合反演能够有效提高高阶部分的信噪比，即提高了反演结果的空间分辨率。

为了验证本发明的技术效果，可以根据步骤3中联合反演模型反演局部地表质量变化，与单独利用GRACE或GNSS反演结果进行比较，以反演结果与真实信号残差的标准差以及阶方差RMS作为反演结果精度的评判。

本实施例中，利用GRACE、GNSS独立反演的局部地表质量变化与原始信号残差的标准差分别为37.86mm和36.79mm，依据本发明联合GRACE和GNSS反演的局部地表质量变化与原始信号残差的标准差为29.40mm，明显小于两类数据单独反演结果的残差标准差，表明联合反演结果的精度更优。从阶方差RMS及其信噪比分析可以看出，GRACE和GNSS联合反演结果在高阶部分优于单独使用GRACE或GNSS的反演结果，表明联合反演结果的空间分辨率更高。

综上所述，本发明提供了一种联合GRACE星间重力位差和GNSS垂直位移反演精细局部地表质量变化的方法。主要特征表现在，利用GRACE和GNSS在反演局部地表质量变化各自的优势，实现了优势互补，两类数据联合反演可以提高局部地表质量变化估计结果的精细度和可靠性。最终通过联合反演方法在南美洲进行数值模拟测试，验证了本发明的正确性和有效性。

具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。

在一些可能的实施例中，提供一种联合GRACE和GNSS的局部地表质量变化反演系统，包括以下模块，

第一模块，用于获取由局部地表质量变化引起的GRACE星间重力位差和GNSS垂直位移观测值；

第二模块，用于根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程，得到相应法方程；根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程，得到相应法方程；

第三模块，用于利用第二模块所得结果根据最小二乘联合平差反演局部地表质量变化，包括根据法方程结合先验信息方程形成GRACE和GNSS联合反演模型，然后给定两类观测值噪声方差以及正则化参数初始值，利用方差分量估计通过迭代计算确定最优权比，得到联合反演的局部地表质量变化。

在一些可能的实施例中，提供一种联合GRACE和GNSS的局部地表质量变化反演系统，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种联合GRACE和GNSS的局部地表质量变化反演方法。

在一些可能的实施例中，提供一种联合GRACE和GNSS的局部地表质量变化反演系统，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种联合GRACE和GNSS的局部地表质量变化反演方法。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。