联合GNSS三维形变的区域陆地水储量变化反演方法及系统转让专利
申请号 : CN202111084343.7
文献号 : CN113899301B
文献日 : 2022-07-15
发明人 : 钟波 , 李贤炮
申请人 : 武汉大学
摘要 :
本发明提供一种联合GNSS三维形变的区域陆地水储量变化反演方法及系统，包括根据质量负荷的格林函数理论，建立GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系；根据质量负荷格林函数理论建立区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程，并得到相应法方程；根据最小二乘联合平差反演区域陆地水储量变化，包括根据法方程并结合先验信息方程形成GNSS垂直和水平形变的联合反演模型，然后给定观测值噪声方差以及正则化参数初始值，利用方差分量估计方法迭代计算各类数据的最优权比，得到最优联合反演的区域陆地水储量变化。本发明提高了区域陆地水储量变化反演结果的精度和可靠性。
权利要求 :
1.一种联合GNSS三维形变的区域陆地水储量变化反演方法，其特征在于：包括如下步骤，步骤1，根据质量负荷的格林函数理论，建立GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系；
步骤2，根据质量负荷的格林函数理论，建立区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程，然后根据GNSS北向、东向和垂直形变的观测方程形成联合反演法方程组；
步骤3，利用步骤2所得结果根据最小二乘联合平差反演区域陆地水储量变化，包括根据法方程并结合先验约束方程形成GNSS三维形变联合反演模型，然后给定观测值噪声方差以及正则化参数初始值，利用方差分量估计方法迭代计算各类数据的最优权比，得到联合反演的区域陆地水储量变化。
2.如权利要求1所述的联合GNSS三维形变的区域陆地水储量变化反演方法，其特征在于：步骤1中，GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系可以表示为；
其中，n(θ)、e(θ)和u(θ)分别为GNSS北向、东向和垂直形变；ΔM为水文质量负荷，MR为地球质量；R为地球平均半径；ln和hn为负荷Love数；θ为格网点与GNSS测站之间的角距；pn为勒让德多项式；α为水平形变与北向之间的夹角。
3.如权利要求1所述的联合GNSS三维形变的区域陆地水储量变化反演方法，其特征在于：步骤2中区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程的表达式如下；
其中，yN、yE和yU分别对应GNSS北向、东向和垂直形变观测值，AN、AE和AU为三维形变对应的格林函数设计矩阵，x是待估的陆地水储量变化参数，eN、eE和eU分别对应北向、东向和垂直形变观测值残差向量， 和 是观测值误差方差，IN、IE和IU是和GNSS北向、东向和垂直形变观测值yN、yE和yU相关的单位矩阵。
4.如权利要求3所述的联合GNSS三维形变的区域陆地水储量变化反演方法，其特征在于：步骤2中联合反演法方程组的具体表达式如下，其中， 是AN的转置，PN是北向形变观测值yN对应的权阵，
是AE的转置，PE是东向形变观测值yE对应的权阵，
是AU的转置，PU是垂直形变观测值yU对应的权阵，
5.如权利要求1所述的联合GNSS三维形变的区域陆地水储量变化反演方法，其特征在于：步骤3中，先验约束方程为Lx＝0，利用离散拉普拉斯矩阵L对反演结果进行正则化约束；
其中，L为拉普拉斯矩阵，具体由拉普拉斯模板 L2＝[1 ‑2 1]和组成，x是待估的陆地水储量变化参数，0是和观测值长度一致的零向量。
6.如权利要求4所述的联合GNSS三维形变的区域陆地水储量变化反演方法，其特征在于：步骤3中根据法方程并结合先验约束方程形成GNSS三维形变联合反演模型，即形成联合
2 2 2 2
反演目标函数：||ANx‑yN||+||AEx‑yE|| +||AUx‑yU|| +λ||Lx||＝min，其中λ为正则化参数，L为拉普拉斯矩阵；
则待求的区域陆地水储量变化参数估计值的联合反演模型可以表示为：
其中，PN、PE和PU分别为北向形变观测值yN、东向形变观测值yE和垂向形变观测值yU对应的权阵；
给定三类观测值噪声方差 以及正则化参数α一个初始值，通过方差分
量估计迭代定权求解最优的区域陆地水储量变化结果。
7.联合GNSS三维形变的区域陆地水储量变化反演系统，其特征在于：包括以下模块，第一模块，用于根据质量负荷的格林函数理论，建立GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系；
第二模块，用于根据质量负荷的格林函数理论，建立区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程，然后根据GNSS北向、东向和垂直形变观测方程形成联合反演法方程组；
第三模块，用于利用第二模块所得结果根据最小二乘联合平差反演区域陆地水储量变化，包括根据法方程并结合先验约束方程形成GNSS三维形变联合反演模型，然后给定观测值噪声方差以及正则化参数初始值，利用方差分量估计方法迭代计算各类数据的最优权比，得到联合反演的区域陆地水储量变化。
8.一种计算机设备，包括处理器和存储器，所述存储器存储有计算机程序，其特征在于：所述处理器执行所述计算机程序时实现权利要求1‑6中任一项方法的步骤。
9.一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，其特征在于：所述计算机程序被处理器执行时实现权利要求1‑6所述的任一项方法的步骤。
说明书 :
联合GNSS三维形变的区域陆地水储量变化反演方法及系统
技术领域
[0001] 本发明涉及利用空间大地测量观测数据反演区域陆地水储量变化的方法，特别涉及一种联合GNSS三维形变的区域陆地水储量变化反演方法及系统。
背景技术
[0002] 陆地水储量(Terrestrial Water Storage,TWS)是指存储在地表以及地下的全部水分，包括地表水、地下水、冰川、积雪及土壤水等，是全球水循环的重要组成部分。精确测定陆地水储量变化的迁移与重新分布对自然环境变化及人类活动具有重要影响。根据负荷格林函数理论，地表水负荷变化会导致地壳发生形变。因此，通过观测高精度的地表形变场及时空变化可以推演和监测地球系统的物质迁移和交换过程，并且获取地表形变场的精细程度越高，其包含的地球系统物质迁移信息则越丰富，这对大地测量学、水文学、冰川学以[1‑3]及全球环境变化等研究都具有十分重要的物理意义 。
[0003] 自Argus et al.(2014)首次提出利用GNSS垂直形变反演区域地表质量变化以来，GNSS 反演的地表质量变化被认为是GRACE/GFO反演结果的有益补充。当前大多数学者主要是利用GRACE/GFO来反演全球和区域尺度上的陆地水储量变化，但受限于其卫星轨道高度、有效载荷以及各类误差源的影响，其反演结果仍然存在诸多局限性。陆地水储量变化在导致地球外部重力场变化的同时也会导致地壳发生形变，这种形变能够被GNSS以毫米级精度[1,5]连续地观测到，因此利用GNSS地表形变资料可独立地反演陆地水储量变化 。目前，数据公开的全球GNSS连续站已达数千个，很多国家还建立了密集的GNSS连续运行台网(如中国[6‑7]
“陆态网络”)对区域地表形变和地质构造运动进行监测 。近年来，随着GNSS观测数据的积累以及数据处理方法的改善，这些高精度、丰富的GNSS数据为研究区域陆地水储量变化提供了重要的数据支撑。
“陆态网络”)对区域地表形变和地质构造运动进行监测 。近年来，随着GNSS观测数据的积累以及数据处理方法的改善，这些高精度、丰富的GNSS数据为研究区域陆地水储量变化提供了重要的数据支撑。
[0004] 现有的研究成果侧重于单独利用GNSS垂直形变反演区域陆地水储量变化，并将其[8]与 GRACE/GFO以及水文模型反演结果进行对比分析 。Wahr et al.(2013)的研究发现，将GNSS 水平形变和垂直形变结合起来能够更好地确定质量负荷的空间分布。在此基础上，本发明创新性地提出联合GNSS三维形变反演区域陆地水储量变化，以进一步提高区域陆地水储量变化反演结果的可靠性。
[0005] 相关参考文献如下：
[0006] [1]Blewitt G.,D.Lavallée,P.Clarke,K.Nurutdinov.A new global mode of earth deformation: seasonal cycle detected.Science,2001,294(5550):2342‑2345.[0007] [2]Tapley,B.D.；Bettadpur,S.；Ries,J.C.；Thompson,P.F.；Watkins,M.M.GRACE Measurements of Mass Variability in the Earth System.Science 2004,305,503–505.
[0008] [3]许厚泽,陆洋,钟敏,郑伟,张子占.卫星重力测量及其在地球物理环境变化监测中的应用.中国科学:地球科学,2012,42(6):843‑853.
[0009] [4]Argus,D.F.；Fu,Y.；Landerer,F.Seasonal variation in total water storage in California inferred from GPS observations of vertical land motion.Geophys.Res.Lett.2014,41,1971–1980.
[0010] [5]Wu,X.Large‑scale global surface mass variations inferred from GPS measurements of load‑induced deformation[J].Geophysical Research Letters,2003,30(14):253‑266.
[0011] [6]李强,游新兆,杨少敏,等.中国大陆构造变形高精度大密度GPS监测—现今速度场[J]. 中国科学:地球科学,2012,42(005):629‑632.
[0012] [7]姜卫平.卫星导航定位基准站网的发展现状、机遇与挑战[J].测绘学报, 2017(10):181‑190.
[0013] [8]Zhong B,Li X,Chen J,et al.Surface Mass Variations from GPS and GRACE/GFO:A Case Study in Southwest China[J].Remote Sensing,2020,12(11):1835.[0014] [9]Wahr,J.；Khan,S.A.；Van Dam,T.；Liu,L.；Van Angelen,J.H.；Broeke,M.R.V.D.；Meertens, C.The use of GPS horizontals for loading studies,with applications to northern California and southeastGreenland.J.Geophys.Res.Solid Earth 2013,118,1795–1806.
发明内容
[0015] 针对现有单独利用GNSS垂直形变反演区域陆地水储量变化方法的不足，本发明提供了一种联合GNSS三维形变反演区域陆地水储量变化的方法及系统。
[0016] 本发明采用的技术方案为一种联合GNSS三维形变反演区域陆地水储量变化的反演方法，包括如下步骤，
[0017] 步骤1，根据质量负荷的格林函数理论，建立GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系；
[0018] 步骤2，根据质量负荷的格林函数理论，建立区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程，然后根据GNSS垂直和北向、东向观测方程形成联合反演法方程组；
[0019] 步骤3，利用步骤2所得结果根据最小二乘联合平差反演区域陆地水储量变化，包括根据法方程并结合先验约束方程形成GNSS三维形变联合反演模型，然后给定观测值噪声方差以及正则化参数初始值，利用方差分量估计方法迭代计算各类数据的最优权比，得到联合反演的区域陆地水储量变化。
[0020] 进一步的，步骤1中，GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系可以表示为；
[0021]
[0022] 其中，n(θ)、e(θ)和u(θ)分别为GNSS北向、东向和垂直形变；ΔM为水文质量负荷，MR为地球质量；R为地球平均半径；ln和hn为负荷Love数；θ为格网点与GNSS测站之间的角距；pn为勒让德多项式；α为水平形变与北向之间的夹角。
[0023] 进一步的，步骤2中区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程的表达式如下；
[0024]
[0025] 其中，yN、yE和yU分别对应GNSS北向、东向和垂直形变观测值，AN、AE和AU为三维形变对应的格林函数设计矩阵，x是待估的陆地水储量变化参数，eN、eE和eU分别对应北向、东向和垂直形变观测值残差向量， 和 是观测值误差方差，IN、IE和IU是和GNSS北向、东向和垂直形变观测值yN、yE和yU相关的单位矩阵。
[0026] 进一步的，步骤2中联合反演法方程组的具体表达式如下，
[0027]
[0028] 其中， 是AN的转置，PN是北向形变观测值yN对应的权阵，是AE的转置，PE是东向形变观测值yE对应的权阵，
是AU的转置，PU是垂直形变观测值yU对应的权阵，
是AU的转置，PU是垂直形变观测值yU对应的权阵，
[0029] 进一步的，步骤3中，先验约束方程为Lx＝0，利用离散拉普拉斯矩阵L对反演结果进行正则化约束；
[0030] 其中，L为拉普拉斯矩阵，具体由拉普拉斯模板 L2＝[1 ‑2 1]和组成，x是待估的陆地水储量变化参数，0是和观测值长度一致的零向量；
[0031] 进一步的，步骤3中根据法方程并结合先验约束方程形成GNSS三维形变联合反演2 2 2 2
模型，即形成联合反演目标函数：||ANx‑yN||+||AEx‑yE||+||AUx‑yU||+λ||Lx||＝min，其中λ为正则化参数，L为拉普拉斯矩阵；
模型，即形成联合反演目标函数：||ANx‑yN||+||AEx‑yE||+||AUx‑yU||+λ||Lx||＝min，其中λ为正则化参数，L为拉普拉斯矩阵；
[0032] 则待求的区域陆地水储量变化参数估计值的联合反演模型可以表示为：
[0033]
[0034] 其中，PN、PE和PU分别为北向形变观测值yN、东向形变观测值yE和垂向形变观测值yU对应的权阵；
[0035] 给定三类观测值噪声方差 以及正则化参数α一个初始值，通过方差分量估计迭代定权求解最优的区域陆地水储量变化结果。
[0036] 另一方面，本发明还提供一种联合GNSS三维形变的区域陆地水储量变化反演系统，用于实现如上所述的一种联合GNSS三维形变的区域陆地水储量变化反演方法，包括以下模块，
[0037] 第一模块，用于根据质量负荷的格林函数理论，建立GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系；
[0038] 第二模块，用于根据质量负荷的格林函数理论，建立区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程，然后根据GNSS垂直和北向、东向观测方程形成联合反演法方程组；
[0039] 第三模块，用于利用第二模块所得结果根据最小二乘联合平差反演区域陆地水储量变化，包括根据法方程并结合先验约束方程形成GNSS三维形变联合反演模型，然后给定观测值噪声方差以及正则化参数初始值，利用方差分量估计方法迭代计算各类数据的最优权比，得到联合反演的区域陆地水储量变化。
[0040] 本发明还提供一种计算机设备，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种联合GNSS三维形变的区域陆地水储量变化反演方法。
[0041] 本发明还提供一种计算机可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种联合GNSS三维形变的区域陆地水储量变化反演方法。
[0042] 与现有技术相比，本发明具有以下优点和有益效果：
[0043] 1、提高了区域陆地水储量变化反演的可靠性。
[0044] 受限于GNSS测站分布稀疏及分布范围有限等因素，当前单独利用GNSS垂直形变数据反演区域陆地水储量变化存在一定的缺陷与不足。本发明联合GNSS三维形变数据反演区域陆地水储量变化的方法，可以充分地利用GNSS三维形变观测值，以对陆地水储量变化的空间分布进行更好的约束，进而提高了反演结果的可靠性。
[0045] 2、方案的理论依据充分，实现稳定有效。
[0046] 通过分别对GNSS北向、东向和垂直形变观测值反演结果以及GNSS三维形变联合反演结果与原始信号的残差标准差进行对比可知，联合反演结果与原始信号残差更小，反映出联合反演结果可靠性更高。本发明充分利用GNSS三维形变观测数据，因而联合反演方法对提高区域陆地水储量变化的可靠性是有效的。
附图说明
[0047] 图1为本发明实施例的方法流程图。
具体实施方式
[0048] 以下结合附图和实施例对本发明技术方案进行具体描述。
[0049] 本发明充分利用GNSS三维形变反演区域陆地水储量变化的优缺点，实现优势互补，通过联合GNSS三维形变观测数据反演可靠的区域陆地水储量变化。由此本发明提出了一种联合GNSS三维形变数据反演可靠的区域陆地水储量变化的方法。
[0050] 下面以中国大陆地区作为实验区域，选取由GLDAS计算得到的2005年9月的水储量变化作为原始信号，通过闭合数值模拟并结合图1说明本发明实施例提供的一种联合GNSS三维形变数据的区域陆地水储量变化反演方法。实施例中，为了方便计算，研究区域被划分为1 °× 1 °的均匀格网，包括海洋区域在内共有3240个格网。研究区所选GNSS测站数为272个，利用GLDAS水文模型计算的水储量变化原始信号和GNSS测站经纬度，根据研究区格网经纬度由
和 单位为m，可以建立GNSS北向、东向和垂直形变与区域
陆地水储量变化之间的数学关系。GNSS北向、东向和垂直形变中分别加入噪声标准差分别为0.0001m、0.0001m以及0.0005m的高斯随机白噪声。具体流程步骤实现如下：
和 单位为m，可以建立GNSS北向、东向和垂直形变与区域
陆地水储量变化之间的数学关系。GNSS北向、东向和垂直形变中分别加入噪声标准差分别为0.0001m、0.0001m以及0.0005m的高斯随机白噪声。具体流程步骤实现如下：
[0051] 步骤1，根据质量负荷的格林函数理论，GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系可以表示为：
[0052]
[0053] 其中，n(θ)、e(θ)和u(θ)分别为GNSS北向、东向和垂直形变；ΔM为水文质量负荷，MR为地球质量；R为地球平均半径；ln和hn为负荷Love数；θ为格网点与GNSS测站之间的角距；pn为勒让德多项式；α为水平形变与北向之间的夹角。
[0054] 步骤2，根据质量负荷的格林函数理论建立区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程；根据GNSS垂直和水平形变(北向、东向)观测方程形成联合反演法方程组。
[0055] GNSS北向、东向和垂直形变的观测方程为：
[0056]
[0057] 其中，yN、yE和yU分别对应GNSS北向、东向和垂直形变观测值，AN、AE和AU为三维形变对应的格林函数设计矩阵，x是待估的陆地水储量变化参数，eN、eE和eU分别对应北向、东向和垂直形变观测值残差向量， 和 是观测值误差方差，IN、IE和IU是和GNSS北向、东向和垂直形变观测值yN、yE和yU相关的单位矩阵。
[0058] 根据GNSS北向、东向和垂直形变的观测方程便可得到其法方程组：
[0059]
[0060] 其中， 是AN的转置，PN是北向形变观测值yN对应的权阵，是AE的转置，PE是东向形变观测值yE对应的权阵，
是AU的转置，PU是垂直形变观测值yU对应的权阵，
是AU的转置，PU是垂直形变观测值yU对应的权阵，
[0061] 为了验证联合反演相对于单独使用GNSS北向、东向和垂直形变反演区域陆地水储量变化的结果更优，根据上述单个观测方程并结合先验约束方程便可以得到单独使用GNSS北向、东向或垂直形变反演区域陆地水储量变化的反演方程：
[0062]
[0063] 其中，σ2为噪声方差，λ为正则化参数，A为设计矩阵，y表示观测值，min表示最小 (具体可理解为当观测值与模型以及约束矩阵的残差二范数最小时则可以求出最优的待估的区域水储量变化参数x)。
[0064] 待求的水储量变化 可以表示为：
[0065]
[0066] 这里P为观测值的权阵，当给定噪声方差σ2和正则化参数λ一个初值，便可以通过方差分量估计迭代计算出区域陆地水储量变化。
[0067] 实施例中，GNSS北向形变的观测方程为： GNSS东向形变的观测方程为： GNSS垂直形变的观测方程为：
先验约束方程为Lx＝0。分别利用含误差的GNSS三维形变模
拟观测数据，根据单一类型数据反演模型 采用方差分量估计
可以迭代计算出区域陆地水储量变化。可以分别得到GNSS北向、东向和垂直形变的反演结果，以及北向、东向和垂直形变反演结果与原始信号的残差图。GNSS垂直形变对水文负荷信号的变化更为敏感，但是其观测精度较低。GNSS北向和东向形变观测精度较高，但其对水文质量负荷变化的敏感性较低。因此，可以综合利用GNSS三维形变数据实现优势互补，进而反演出最佳的区域陆地水储量变化。此外，在模拟研究区域(中国大陆地区)，GNSS北向、东向和垂直形变独立反演结果与原始信号残差的标准差分别为24.93mm、30.32mm以及
26.11mm。
先验约束方程为Lx＝0。分别利用含误差的GNSS三维形变模
拟观测数据，根据单一类型数据反演模型 采用方差分量估计
可以迭代计算出区域陆地水储量变化。可以分别得到GNSS北向、东向和垂直形变的反演结果，以及北向、东向和垂直形变反演结果与原始信号的残差图。GNSS垂直形变对水文负荷信号的变化更为敏感，但是其观测精度较低。GNSS北向和东向形变观测精度较高，但其对水文质量负荷变化的敏感性较低。因此，可以综合利用GNSS三维形变数据实现优势互补，进而反演出最佳的区域陆地水储量变化。此外，在模拟研究区域(中国大陆地区)，GNSS北向、东向和垂直形变独立反演结果与原始信号残差的标准差分别为24.93mm、30.32mm以及
26.11mm。
[0068] 步骤3，利用步骤2所得结果根据最小二乘联合平差反演区域陆地水储量变化，包括根据法方程并结合先验约束方程形成GNSS三维形变联合反演模型，然后给定观测值噪声方差以及正则化参数初始值，利用方差分量估计方法迭代计算各类数据的最优权比，得到联合反演的区域陆地水储量变化。
[0069] 由于GNSS三维形变数据对水文质量负荷变化的敏感性差异以及观测精度的不同，要实现这两类观测值的最佳融合，需要充分考虑其最优融合的模式及准则，本发明通过方差分量估计确定两类数据的最优权比。
[0070] 联合GNSS三维形变反演区域陆地水储量变化属于离散病态问题，需要利用恰当的正则化方法对病态法方程组进行约束求解，常用的包括Tikhonov正则化、截断奇异值分解(TSVD) 以及岭估计等方法都受限于如何快速有效地获取最优正则化参数，本发明以拉普拉斯矩阵作为先验空间约束矩阵，通过最小二乘迭代方法由观测数据本身自适应地、有效地获取正则化参数。
[0071] 本发明进一步提出，联合反演需要确定联合反演模型和最优权比，还需包含以下子步骤：
[0072] 3.1根据法方程并结合先验信息方程形成GNSS三维形变联合反演模型，即形成联2 2 2 2
合反演目标函数：||ANx‑yN||+||AEx‑yE||+||AUx‑yU||+λ||Lx||＝min，其中λ为正则化参数。
合反演目标函数：||ANx‑yN||+||AEx‑yE||+||AUx‑yU||+λ||Lx||＝min，其中λ为正则化参数。
[0073] 则待求的区域陆地水储量变化参数估计值的联合反演模型可以表示为：
[0074] 3.2给定观测值噪声方差以及正则化参数初始值，利用方差分量估计方法通过迭代计算确定最优权比，可以最优联合反演的区域陆地水储量变化。
[0075] 给定 以及正则化参数λ一个初始值，通过方差分量估计迭代定权便可以估计最优的区域陆地水储量变化结果。
[0076] 本实施例中，分别利用GNSS北向、东向、垂直形变以及联合三维形变反演结果与原始信号的残差标准差进行统计分析便可以评估联合反演结果可靠性。
[0077] 实施例中，联合带有误差的GNSS三维形变模拟数据反演中国大陆地区陆地水储量变化。
[0078] 具体包含以下子步骤：
[0079] 3.1根据步骤2中GNSS东向、北向和垂直形变的观测方程及先验约束方程则可以得到联合反演模型：
[0080]
[0081] 由于利用GNSS三维形变反演区域陆地水储量变化属于离散病态问题，利用离散拉普拉斯矩阵L对反演结果进行正则化约束，其约束方程为：
[0082] Lx＝0，
[0083] 其中，L为拉普拉斯矩阵，具体由拉普拉斯模板 L2＝[1 ‑2 1]和组成，x是待估的陆地水储量变化参数，0是和观测值长度一致的零向量；
[0084] 3.2以模拟的GNSS三维形变的噪声标准差作为方差初始值并给定一个正则化参数，通过方差分量估计迭代定权估计最优的区域陆地水储量变化结果。得到GNSS三维形变联合反演结果，由于GNSS垂直形变对水文负荷信号的变化更为敏感，但是其观测精度较低；GNSS 北向和东向形变观测精度较高，但其对水文质量负荷变化的敏感性较低；相对于单独利用 GNSS北向、东向或垂直形变反演陆地水储量变化而言，联合反演可以综合利用三维形变信息对陆地水储量变化空间分布进行更好的约束，进而得到更加可靠的区域陆地水储量变化反演结果。可得到单独利用GNSS北向、东向或垂直形变反演与原始信号的残差图，其中反演结果与原始信号残差标准差为24.93mm、30.32mm和26.11mm。可得到联合反演结果与原始信号的残差图，其中反演结果与原始信号残差标准差为18.52mm。根据结果可知，联合反演结果与原始信号残差的STD远小于单独利用一维形变数据反演结果与真实信号残差的STD，联合反演结果能够充分利用GNSS形变观测数据，进而反演得到更加可靠的区域水储量变化反演结果。
[0085] 为了验证本发明的技术效果，可以根据步骤3中联合反演模型解算区域陆地水储量变化，并与单独利用GNSS北向、东向和垂直形变反演结果进行比较，以反演结果与真实信号残差的标准差作为反演结果可靠性的评判。
[0086] 本实施例中，单独利用GNSS北向、东向和垂直形变反演的区域陆地水储量变化与原始信号残差的标准差分别为24.93mm、30.32mm和26.11mm，依据本发明联合GNSS反演的区域陆地水储量变化与原始信号残差的标准差为18.52mm，明显小于两类数据单独反演结果的残差标准差，表明联合反演结果的精度和可靠性更优。
[0087] 综上所述，本发明提供了一种联合GNSS三维形变反演区域陆地水储量变化的方法。主要特征表现在，综合利用GNSS北向、东向和垂直形变在反演区域陆地水储量变化的各自优势，有效实现了优势互补，并且联合反演可以提高区域陆地水储量变化估计结果的精度和可靠性。最终通过联合反演方法在中国大陆地区进行数值模拟测试，验证了本发明的正确性和有效性。
[0088] 具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。
[0089] 在一些可能的实施例中，提供一种联合GNSS三维形变反演区域陆地水储量变化反演系统，包括以下模块，
[0090] 第一模块，用于根据质量负荷的格林函数理论，建立GNSS北向、东向和垂直形变与区域陆地水储量变化之间的关系；
[0091] 第二模块，用于根据质量负荷的格林函数理论，建立区域陆地水储量变化与GNSS北向、东向和垂直形变之间的观测方程，然后根据GNSS垂直和北向、东向观测方程形成联合反演法方程组；
[0092] 第三模块，用于利用第二模块所得结果根据最小二乘联合平差反演区域陆地水储量变化，包括根据法方程并结合先验约束方程形成GNSS三维形变联合反演模型，然后给定观测值噪声方差以及正则化参数初始值，利用方差分量估计方法迭代计算各类数据的最优权比，得到联合反演的区域陆地水储量变化。
[0093] 在一些可能的实施例中，提供一种计算机设备，包括处理器和存储器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现如上所述的一种联合三维形变反演区域陆地水储量变化反演方法。
[0094] 在一些可能的实施例中，提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如上所述的一种联合三维形变反演区域陆地水储量变化反演方法。
[0095] 本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。