CN112989589B - 联合grace和gnss的局部地表质量变化反演方法及系统 - Google Patents

联合grace和gnss的局部地表质量变化反演方法及系统 Download PDF

Info

Publication number
CN112989589B
CN112989589B CN202110244382.2A CN202110244382A CN112989589B CN 112989589 B CN112989589 B CN 112989589B CN 202110244382 A CN202110244382 A CN 202110244382A CN 112989589 B CN112989589 B CN 112989589B
Authority
CN
China
Prior art keywords
grace
gnss
surface quality
earth surface
local
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110244382.2A
Other languages
English (en)
Other versions
CN112989589A (zh
Inventor
钟波
李贤炮
谭江涛
刘滔
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202110244382.2A priority Critical patent/CN112989589B/zh
Publication of CN112989589A publication Critical patent/CN112989589A/zh
Application granted granted Critical
Publication of CN112989589B publication Critical patent/CN112989589B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10—Noise analysis or noise optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Computing Systems (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

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

Description

联合GRACE和GNSS的局部地表质量变化反演方法及系统
技术领域
本发明涉及联合空间大地测量观测数据反演局部地表质量变化的方法，特别涉及一种联合GRACE星间重力位差和GNSS垂直形变数据反演精细局部地表质量的方法及系统。
背景技术
地表质量变化主要包括陆地水储量变化、极地冰盖和山岳冰川融化、大气压变化、海洋质量变化以及和固体地球物理现象相关的其他质量变化现象。这些全球或局部的地表质量迁移与重新分布，不仅会引起地球重力场的改变，还会导致固体地球表面产生形变。因此，通过观测地球重力场和地球形变场的时空变化可以推演和监测地球系统的物质迁移和交换过程，并且获取时变重力场与形变场的精细程度越高，其包含的地球系统物质迁移信息则越丰富，这对大地测量学、固体地球物理学、地球动力学、海洋学、冰川学以及全球环境变化等研究都具有十分重要的物理意义[1-3]。
近年来，卫星重力测量技术的实现，特别是美德联合实施的GRACE计划，为准实时获取全球高精度时变重力场进而反演地表质量变化提供了新途径。不同于传统的重力探测技术，GRACE能够获取全球高精度、均匀覆盖的月时变重力场，可在季节时间尺度上以1cm等效水高精度、400km空间分辨率提供地表质量变化信息[2]，实现了定量揭示全球环境变化导致地表质量迁移的里程碑式突破。尽管GRACE卫星重力测量已经证明了以前所未有的精度测量与陆地水文循环、冰川融化、海平面变化、地震和人类活动相关的大规模质量变化的巨大潜力，但其仍然存在很大的局限性：①GRACE球谐系数的截断以及需要进行空间滤波和平滑处理以抑制高频和南北条带噪声的过程会导致信号泄漏，进而降低了GRACE反演地表质量变化的空间分辨率(～300–500km)[4]；②GRACE低阶球谐系数的不确定性相对较大，特别是二阶带谐项C20，需要利用诸如卫星激光测距(SLR)的其他大地测量反演结果来代替[5]；③2016年8月之后，GRACE观测数据质量下降，影响了GRACE反演地表质量变化的准确性；④GRACE和其继任者GRACE Follow-On(GFO)之间大约有一年的数据间隔，GRACE/GFO数据的不连续性影响了其科学解释和应用。
由于地球质量重新分布导致的地表形变能够被GNSS以毫米级精度连续地观测，因此可以利用GNSS负荷形变资料独立地反演地表质量变化[1,6]。自从Blewitt等(2001)提出地表负荷反演方法以来，利用GNSS地表形变数据反演的地表质量变化一直是全球和区域尺度上GRACE/GFO观测结果的补充[1,7]。相较于GRACE/GFO，GNSS有如下优势：①GNSS网络在全球范围内分布广泛，并且可以以mm级的精度和较高的时空分辨率近乎实时地获取观测数据；②GNSS地表形变(特别是垂直位移)，对局部(～10km)和区域质量负荷的变化更为敏感，通过密集的GNSS观测网络反演的局部地表质量变化的空间分辨率可以达到～50km[8-9]；③利用密集的GNSS观测资料还能以更高的时空分辨率反演局部地表质量变化，有助于研究局部短时间尺度的地球动力学过程。但受限于自然条件和地区经济发展的影响，GNSS测站分布数量仍然相对有限，在广阔的海洋区域和自然环境恶劣的地区难以布设GNSS站点，或者在经济落后的国家或地区难以支撑高昂的站点布设与管理费用等，这些因素都制约了GNSS技术的进一步应用。
目前国内外学者的研究成果侧重于单独使用GRACE或GNSS数据反演局部地表质量变化，或对两类数据的反演结果进行一致性和相关性分析。虽然Adusumilli等(2019)在反演美国周边地区水储量变化时利用GRACE和GNSS观测资料进行联合反演，但其实质是使用GRACE反演结果对GNSS反演结果进行空间约束，并没有从观测数据层面进行联合反演[9]。所以，至今尚未有从GRACE L1B原始观测数据入手联合GNSS垂直形变反演局部地表质量变化的方法。本发明创新性地提出利用GRACE L1B观测数据通过能量平衡方程计算得到星间重力位差数据[11-12]，将其与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星间重力位差之间的观测方程如下，
Figure BDA0002963546610000041
其中，y1是GRACE星间重力位差观测值，A1是星间重力位差设计矩阵，x是待估的地表质量变化参数，e1是重力位差的观测值残差向量，
Figure BDA0002963546610000042
是观测值误差方差，I1是和重力位差观测值y1相关的单位矩阵；
根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程如下，
Figure BDA0002963546610000043
其中，y2是GNSS垂直位移观测值，A2是格林函数设计矩阵，x是待估的地表质量变化参数，e2是垂直位移的观测值残差向量，
Figure BDA0002963546610000044
是观测值误差方差，I2是和GNSS垂直位移观测值y2相关的单位矩阵。
而且，步骤3中，设增加先验约束方程
x0＝Ixx+e0,e0～(0,Cx)，
其中，x0为先验地球物理信息，x是待估的地表质量变化参数，Ix是和待估参数相对应的设计矩阵，e0是满足期望为0和方差矩阵为Cx的残差向量，其中Cx是根据先验信息计算的空间协方差函数；
采用的联合反演模型为，
Figure BDA0002963546610000045
其中，P1、P2分别为GRACE和GNSS观测值的权阵；
给定噪声方差
Figure BDA0002963546610000046
以及正则化参数α一个初始值，通过方差分量估计迭代定权估计最优的局部地表质量变化结果。
而且，根据步骤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双星中点处的星间重力位差与局部地表质量变化之间的关系可以表示为：
Figure BDA0002963546610000061
其中，
T12为星间重力位差；
G为万有引力常数：6.673×10-11m3kg-1s-2；
kn为负荷Love数；
N为研究区格网点数，j为格网点标号，j＝1,2,…N；
Figure BDA0002963546610000062
和
Figure BDA0002963546610000063
分别表示两颗不同卫星到格网点之间的距离，可以展开为
Figure BDA0002963546610000064
其中R为地球的平均半径，
Figure BDA0002963546610000065
表示
Figure BDA0002963546610000066
或
Figure BDA0002963546610000067
和
Figure BDA0002963546610000068
分别表示两颗卫星到格网点之间的球面角距；Pn()为勒让德多项式，且
Figure BDA0002963546610000069
(r1,θ1,λ1)和(r2,θ2,λ2)为GRACE双星的地心半径、地心余纬和经度；δmj为地表质量负荷。
根据质量负荷的格林函数理论，GNSS垂直位移与局部地表质量变化之间的关系可以表示为：
Figure BDA0002963546610000071
其中，u为GNSS垂直位移；MR为地球质量；hn为负荷Love数；θj为格网点与GNSS测站之间的角距。
实施例中，为了方便计算，研究区域被划分为2°×2°的均匀格网，包括海洋区域在内共有1216个格网。以2005年9月的GLDAS水文模型计算的水储量变化作为原始信号，根据研究区的经纬度范围和GRACE卫星星历数据计算得到GRACE星间重力位差数据，包括海洋区域共有38696个观测值。通过方程
Figure BDA0002963546610000072
单位为m2/s2，可以建立GRACE星间重力位差与局部地表质量变化之间的数学关系。研究区所选GNSS测站数为596个，利用GLDAS水文模型计算的水储量变化原始信号和GNSS测站经纬度，根据研究区格网由
Figure BDA0002963546610000073
单位为m，可以建立GNSS垂直位移与局部地表质量变化之间的数学关系。为了与实际观测数据更加接近，GRACE星间重力位差和GNSS垂直位移数据中分别加入标准差为0.001m2/s2和0.001m的高斯随机白噪声。
步骤2，根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程；根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程；根据GRACE和GNSS观测方程形成联合反演法方程组。
GRACE星间重力位差的观测方程为：
Figure BDA0002963546610000074
其中，y1是GRACE星间重力位差观测值，A1是星间重力位差设计矩阵，x是待估的地表质量变化参数，e1是重力位差的观测值残差向量，
Figure BDA0002963546610000075
是观测值误差方差，I1是和重力位差观测值y1相关的单位矩阵。
根据GRACE星间重力位差的观测方程便可得到其法方程：
Figure BDA0002963546610000081
其中，
Figure BDA0002963546610000082
是A1的转置，P1是和星间重力位差观测值y1相关的权阵；
Figure BDA0002963546610000083
GNSS垂直位移的观测方程为：
Figure BDA0002963546610000084
其中，y2是GNSS垂直位移观测值，A2是格林函数设计矩阵，x是待估的地表质量变化参数，e2是垂直位移的观测值残差向量，
Figure BDA0002963546610000085
是观测值误差方差，I2是和GNSS垂直位移观测值y2相关的单位矩阵。
由于利用GRACE或GNSS数据反演局部地表质量变化属于离散病态问题，可利用先验信息对病态方程进行约束求解。根据GRACE或GNSS的观测方程，增加先验约束方程
x0＝Ixx+e0,e0～(0,Cx)，
其中，x0为先验地球物理信息(如水文模型)，x是待估的地表质量变化参数，Ix是和待估参数相对应的设计矩阵，e0是满足期望为0和方差矩阵为Cx的残差向量，其中Cx是根据先验信息计算的空间协方差函数。
根据GNSS垂直位移的观测方程便可得到其法方程：
Figure BDA0002963546610000086
其中，
Figure BDA0002963546610000087
是A2的转置，P2是和GNSS垂直位移观测值y2相关的权阵；
Figure BDA0002963546610000088
根据GRACE星间重力位差观测值和GNSS垂直位移观测值的法方程便可得到联合反演法方程组：
Figure BDA0002963546610000089
单独使用GRACE或GNSS反演局部地表质量变化的反演方程
||Ax-y||2+α||x-x0||2＝min，
其中，α为正则化参数，A为设计矩阵，y表示观测值，min表示最小(具体可理解为当目标函数最小时则可以求出待估的参数x)。
待求的地表质量变化
Figure BDA0002963546610000091
可以表示为：
Figure BDA0002963546610000092
这里P为权阵，通过给定噪声方差
Figure BDA0002963546610000093
和正则化参数α一个初始值，通过方差分量估计可以迭代计算出局部地表质量变化。
实施例中，GRACE星间重力位差的观测方程为：y1＝A1x+e1,
Figure BDA0002963546610000094
GNSS垂直位移的观测方程为：y2＝A2x+e2,
Figure BDA0002963546610000095
先验约束方程为：x0＝Ixx+e0,e0～(0,Cx)，这里以GLDAS水文模型作为先验约束。分别利用含误差的GRACE或GNSS模拟观测数据，根据单一类型数据反演模型
Figure BDA0002963546610000096
采用方差分量估计可以迭代计算出局部地表质量变化。可得GRACE、GNSS的反演结果，以及GRACE、GNSS独立反演结果与原始信号的残差图。由于GRACE观测数据在整个研究区域内均匀分布，GRACE对整个研究区域内大尺度低频信号恢复得较好，但对于一些局部小尺度高频信号却没有很好地恢复出来。由实验结果可知，在有测站分布的区域，GNSS很好地恢复了原始信号的空间分布，特别是在GNSS测站密集分布的区域，局部高频信号也能够很好地恢复出来。但由于GNSS测站分布的空间范围有限，在研究区北部没有测站分布，则不能反演出该区域地表质量变化分布形态。此外，GRACE、GNSS独立反演结果与原始信号残差的标准差分别为37.86mm和36.79mm。
步骤3，通过步骤2所得GRACE和GNSS观测方程的法方程组根据最小二乘联合平差反演局部地表质量变化。
GRACE为“物理”观测量，GNSS为“几何”观测量，为实现这两类观测值的最佳融合，需要充分顾及它们在时空分辨率、空间覆盖和频谱敏感性等方面的差异，合理考虑其最优融合的模式及准则，本发明通过方差分量估计确定两类数据的最优权比。
联合GRACE和GNSS反演局部地表质量变化属于离散病态问题，需要利用正则化方法对病态法方程组进行约束求解，常用的包括Tikhonov正则化、截断奇异值分解(TSVD)以及岭估计等方法都受限于如何快速有效地获取最优正则化参数，本发明以先验地球物理模型信息(如水文模型)作为约束，通过最小二乘迭代方法根据观测数据本身自适应地、有效地获取正则化参数。
本发明进一步提出，联合反演需要确定联合反演模型和定权，还需包含以下子步骤：
3.1根据法方程并结合先验信息方程形成GRACE和GNSS联合反演模型，即形成联合反演目标函数：||A1x-y1||2+||A2x-y2||2+α||x-x0||2＝min，其中α为正则化参数。
则待求的地表质量变化参数估计值的联合反演模型可以表示为：
Figure BDA0002963546610000101
这里P1、P2分别为GRACE和GNSS观测值的权阵；
3.2给定两类观测值噪声方差以及正则化参数初始值，利用方差分量估计通过迭代计算确定最优权比，可以计算出联合反演的局部地表质量变化。
给定噪声方差
Figure BDA0002963546610000102
以及正则化参数α一个初始值，通过方差分量估计迭代定权便可以估计最优的局部地表质量变化结果。
本实施例中，分别利用GRACE、GNSS以及联合GRACE和GNSS反演结果与原始信号的残差标准差和阶方差RMS及信噪比进行统计分析便可以评估联合反演结果精度。
实施例中，联合GRACE和GNSS含误差模拟数据反演南美洲局部地表质量变化。具体包含以下子步骤：
3.1根据步骤2中GRACE与GNSS的观测测方程及先验约束方程则可以得到联合反演反演模型：
Figure BDA0002963546610000103
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的局部地表质量变化反演方法。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (7)

1.一种联合GRACE和GNSS的局部地表质量变化反演方法，其特征是：包括如下步骤，
步骤1，获取由局部地表质量变化引起的GRACE星间重力位差和GNSS垂直位移观测值；
步骤2，根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程，得到相应法方程；根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程，得到相应法方程；
步骤3，利用步骤2所得结果根据最小二乘联合平差反演局部地表质量变化，包括根据法方程结合先验信息方程形成GRACE和GNSS联合反演模型，然后给定两类观测值噪声方差以及正则化参数初始值，利用方差分量估计通过迭代计算确定最优权比，得到联合反演的局部地表质量变化；
步骤2中，
根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程如下，
Figure FDA0003659308390000011
其中，y1是GRACE星间重力位差观测值，A1是星间重力位差设计矩阵，x是待估的地表质量变化参数，e1是重力位差的观测值残差向量，
Figure FDA0003659308390000012
是观测值误差方差，I1是和重力位差观测值y1相关的单位矩阵；
根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程如下，
Figure FDA0003659308390000013
其中，y2是GNSS垂直位移观测值，A2是格林函数设计矩阵，x是待估的地表质量变化参数，e2是垂直位移的观测值残差向量，
Figure FDA0003659308390000014
是观测值误差方差，I2是和GNSS垂直位移观测值y2相关的单位矩阵。
2.如权利要求1所述的联合GRACE和GNSS的局部地表质量变化反演方法，其特征是：步骤3中，设增加先验约束方程
x0＝Ixx+e0,e0～(0,Cx)，
其中，x0为先验地球物理信息，x是待估的地表质量变化参数，Ix是和待估参数相对应的设计矩阵，e0是满足期望为0和方差矩阵为Cx的残差向量，其中Cx是根据先验信息计算的空间协方差函数；
采用的联合反演模型为，
Figure FDA0003659308390000021
其中，P1、P2分别为GRACE和GNSS观测值的权阵；
给定噪声方差
Figure FDA0003659308390000022
以及正则化参数α一个初始值，通过方差分量估计迭代定权估计最优的局部地表质量变化结果。
3.如权利要求1或2所述的联合GRACE和GNSS的局部地表质量变化反演方法，其特征是：根据步骤3中联合反演模型反演局部地表质量变化，与单独利用GRACE或GNSS反演结果进行比较，以反演结果与真实信号残差的标准差以及阶方差RMS作为反演结果精度的评判。
4.一种联合GRACE和GNSS的局部地表质量变化反演系统，其特征在于：用于实现如权利要求1-3任一项所述的一种联合GRACE和GNSS的局部地表质量变化反演方法。
5.根据权利要求4所述联合GRACE和GNSS的局部地表质量变化反演系统，其特征在于：包括以下模块，
第一模块，用于获取由局部地表质量变化引起的GRACE星间重力位差和GNSS垂直位移观测值；
第二模块，用于根据牛顿万有引力定律建立局部地表质量变化与GRACE星间重力位差之间的观测方程，得到相应法方程；根据质量负荷的格林函数理论建立局部地表质量变化与GNSS垂直位移之间的观测方程，得到相应法方程；
第三模块，用于利用第二模块所得结果根据最小二乘联合平差反演局部地表质量变化，包括根据法方程结合先验信息方程形成GRACE和GNSS联合反演模型，然后给定两类观测值噪声方差以及正则化参数初始值，利用方差分量估计通过迭代计算确定最优权比，得到联合反演的局部地表质量变化。
6.根据权利要求4所述联合GRACE和GNSS的局部地表质量变化反演系统，其特征在于：包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如权利要求1-3任一项所述的一种联合GRACE和GNSS的局部地表质量变化反演方法。
7.根据权利要求4所述联合GRACE和GNSS的局部地表质量变化反演系统，其特征在于：包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如权利要求1-3任一项所述的一种联合GRACE和GNSS的局部地表质量变化反演方法。
CN202110244382.2A 2021-03-05 2021-03-05 联合grace和gnss的局部地表质量变化反演方法及系统 Active CN112989589B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110244382.2A CN112989589B (zh) 2021-03-05 2021-03-05 联合grace和gnss的局部地表质量变化反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110244382.2A CN112989589B (zh) 2021-03-05 2021-03-05 联合grace和gnss的局部地表质量变化反演方法及系统

Publications (2)

Publication Number Publication Date
CN112989589A CN112989589A (zh) 2021-06-18
CN112989589B true CN112989589B (zh) 2022-07-05

Family

ID=76352969

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110244382.2A Active CN112989589B (zh) 2021-03-05 2021-03-05 联合grace和gnss的局部地表质量变化反演方法及系统

Country Status (1)

Country Link
CN (1) CN112989589B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113792450A (zh) * 2021-08-16 2021-12-14 中国空间技术研究院 基于机器学习负载模型提高陆地水储量反演准确性的方法
CN113899301B (zh) * 2021-09-15 2022-07-15 武汉大学 联合gnss三维形变的区域陆地水储量变化反演方法及系统
CN114046774B (zh) * 2022-01-05 2022-04-08 中国测绘科学研究院 综合cors网和多源数据的地面形变连续监测方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103064128B (zh) * 2013-01-06 2015-12-09 中国科学院测量与地球物理研究所 基于星间距离误差模型的地球重力场恢复方法
CN103091721B (zh) * 2013-01-10 2015-03-25 中国科学院测量与地球物理研究所 利用不同轨道倾角卫星联合反演地球重力场的方法
CN103091722B (zh) * 2013-01-22 2015-06-17 中国科学院测量与地球物理研究所 基于载荷误差分析原理的卫星重力反演方法
US20180172860A1 (en) * 2015-08-05 2018-06-21 Halliburton Energy Services, Inc Time-Lapsed Seismic Wavefield Monitoring of Downhole Formations
CN106909722B (zh) * 2017-02-10 2019-07-26 广西壮族自治区气象减灾研究所 一种近地面气温的大面积精准反演方法
CN107102332B (zh) * 2017-05-11 2019-12-06 中南大学 基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法
CN108279442A (zh) * 2018-01-30 2018-07-13 中国国土资源航空物探遥感中心 一种应用于大数据计算的航空重力数据物性层析计算方法

Also Published As

Publication number Publication date
CN112989589A (zh) 2021-06-18

Similar Documents

Publication Publication Date Title
CN112989589B (zh) 联合grace和gnss的局部地表质量变化反演方法及系统
Sun et al. Combining physically based modeling and deep learning for fusing GRACE satellite data: can we learn from mismatch?
Ferreira et al. Space-based observations of crustal deflections for drought characterization in Brazil
Chen Satellite gravimetry and mass transport in the earth system
Watkins et al. Improved methods for observing Earth's time variable mass distribution with GRACE using spherical cap mascons
Cazenave et al. Time-variable gravity from space and present-day mass redistribution in theEarth system
Feng et al. Groundwater storage variations in the North China Plain from GRACE with spatial constraints
CN109313268A (zh) 具有宽巷偏差校正值和窄巷偏差校正值的导航卫星的轨道和低延时时钟的确定
CN101403790A (zh) 单频gps接收机的精密单点定位方法
CN104597471A (zh) 面向时钟同步多天线gnss接收机的定向测姿方法
CN109154670A (zh) 导航卫星宽巷偏差确定系统和方法
CN109154664A (zh) 具有低延时时钟校正值的导航卫星轨道和时钟的确定
Han et al. GPS recovery of daily hydrologic and atmospheric mass variation: A methodology and results from the Australian continent
WO2023197714A1 (zh) 一种适用于动态载体平台的gnss多路径误差削弱方法
Wahr Time-variable gravity from satellites
Karegar et al. A new hybrid method for estimating hydrologically induced vertical deformation from GRACE and a hydrological model: an example from Central North America
HAN et al. Time-variable gravity field determination using Slepian functions and terrestrial measurements: a case study in North China with data from 2011 to 2013
Han et al. Novel Along‐Track Processing of GRACE Follow‐On Laser Ranging Measurements Found Abrupt Water Storage Increase and Land Subsidence During the 2021 March Australian Flooding
Li et al. Inversion of GNSS vertical displacements for terrestrial water storage changes using Slepian basis functions
Xiang et al. Characterizing the seasonal hydrological loading over the asian continent using GPS, GRACE, and hydrological model
Razeghi et al. A joint analysis of GPS displacement and GRACE geopotential data for simultaneous estimation of geocenter motion and gravitational field
Xu et al. Lake level changes determined by Cryosat-2 altimetry data and water-induced loading deformation around Lake Qinghai
Klosko et al. Evaluation and validation of mascon recovery using GRACE KBRR data with independent mass flux estimates in the Mississippi Basin
Yin et al. Comparison of vertical surface deformation estimates derived from space‐based gravimetry, ground‐based GPS, and model‐based hydrologic loading over snow‐dominated watersheds in the United States
Han Seasonal clockwise gyration and tilt of the Australian continent chasing the center of mass of the Earth's system from GPS and GRACE

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant