CN116451440A - 一种顾及比例因子的陆地水储量变化反演方法及智能终端 - Google Patents
一种顾及比例因子的陆地水储量变化反演方法及智能终端 Download PDFInfo
- Publication number
- CN116451440A CN116451440A CN202310283956.6A CN202310283956A CN116451440A CN 116451440 A CN116451440 A CN 116451440A CN 202310283956 A CN202310283956 A CN 202310283956A CN 116451440 A CN116451440 A CN 116451440A
- Authority
- CN
- China
- Prior art keywords
- gnss
- station
- observation
- grace
- displacement
- 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.)
- Pending
Links
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 91
- 238000000034 method Methods 0.000 title claims abstract description 57
- 230000008859 change Effects 0.000 title claims abstract description 35
- 238000006073 displacement reaction Methods 0.000 claims abstract description 147
- 238000007781 pre-processing Methods 0.000 claims abstract description 10
- 239000011159 matrix material Substances 0.000 claims description 32
- 238000004613 tight binding model Methods 0.000 claims description 20
- 238000012880 independent component analysis Methods 0.000 claims description 16
- 238000005315 distribution function Methods 0.000 claims description 11
- 239000000872 buffer Substances 0.000 claims description 10
- 230000008030 elimination Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 claims description 3
- 230000009897 systematic effect Effects 0.000 abstract description 4
- 230000004044 response Effects 0.000 description 14
- 238000010586 diagram Methods 0.000 description 13
- 238000005259 measurement Methods 0.000 description 11
- 230000006870 function Effects 0.000 description 9
- 230000000694 effects Effects 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000011160 research Methods 0.000 description 5
- 230000003595 spectral effect Effects 0.000 description 5
- 238000004590 computer program Methods 0.000 description 4
- 238000009499 grossing Methods 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 238000013459 approach Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000005484 gravity Effects 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 230000035945 sensitivity Effects 0.000 description 3
- 239000013598 vector Substances 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 239000003086 colorant Substances 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 239000003673 groundwater Substances 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000000670 limiting effect Effects 0.000 description 2
- 230000007774 longterm Effects 0.000 description 2
- 230000002265 prevention Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000011084 recovery Methods 0.000 description 2
- 230000002829 reductive effect Effects 0.000 description 2
- 230000001360 synchronised effect Effects 0.000 description 2
- 238000006424 Flood reaction Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000001934 delay Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 239000002352 surface water Substances 0.000 description 1
- 238000012731 temporal analysis Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000000700 time series analysis Methods 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
- 238000010792 warming Methods 0.000 description 1
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
-
- 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)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Testing Or Calibration Of Command Recording Devices (AREA)
Abstract
本发明公开了一种顾及比例因子的陆地水储量变化反演方法，所述方法包括：对GNSS三维坐标时间序列进行预处理；根据处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；根据GRACE‑Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；根据GNSS观测位移和GRACE模型预测位移，得到GNSS真实观测台站位置的比例因子；根据GNSS真实观测台站位置的比例因子，内插得到虚拟台站的比例因子；根据GNSS真实观测台站和虚拟台站的观测值以及所述虚拟台站的比例因子建立联合反演ΔTWS模型，进行反演。本发明依赖于来自虚拟台站的比例因子的约束，从而可以最大程度减少观测手段之间的系统差异，提供更准确的ΔTWS，对稀疏GNSS台站覆盖区域的水文应用和水资源管理有重要意义。
Description
技术领域
本发明涉及水文分析领域，具体涉及一种顾及比例因子的陆地水储量变化反演方法及智能终端。
背景技术
陆地水储量变化(ΔTWS)是指地表水、地下水、土壤水分、雨、雪、植物冠层水等变化在内的总和，ΔTWS时空演变准确估计对研究水文循环和水资源管理至关重要。特别是，在当前全球变暖的背景下，极端水文事件(洪灾和旱灾)发生的风险明显上升，不仅带来了严重的水和粮食安全问题，还可能进一步导致经济风险和区域冲突。监测ΔTWS有助于分析极端气象水文事件发生的频率和强度变化，进而给出可能的应对和预防策略。传统基于点观测的方法(如水文观测站、测井)可以实现对不同的水成分进行独立测量，但由于其时空覆盖较差，无法在全球或区域范围内获取准确的ΔTWS。水文模型被广泛用于估计全球或区域的ΔTWS，以应对和缓解水资源压力。但水文模型并没有纳入对所有水文成分的观察，因此无法区分气候变化和人类活动对于水资源的影响。
目前，随着多种测地手段的发展，密集的全球导航卫星系统(GNSS)网络和重力恢复和气候实验(GRACE)均提供了对于ΔTWS可靠的独立观测。GRACE通过观察地球重力场的变化来监测ΔTWS质量的重新分布，可以显示粗略空间(>300km)上的水文循环过程，为大规模的干旱、洪水监测以及长期地下水变化提供了重要的信息。但受限于GRACE本身粗糙的空间分辨率，在较短的空间距离使用GRACE推断ΔTWS存在较大的不确定性。GNSS可以准确监测局部(0-200km)地表质量的重新分布引起的负载形变，并基于固体地球的弹性载荷理论来推断ΔTWS。随着GNSS定位算法精度和信号处理算法性能的提升，GNSS被越来越多的用于极端水文干旱事件、极端降水事件、冰川质量损失估计以及雪水当量估计等方面的研究。但是，GNSS反演ΔTWS的可靠性依赖于GNSS台站的空间分布和密度，在稀疏GNSS台站覆盖的区域反演性能可能受到影响，这极大地限制了该方法的适用性和准确性。
最近的研究开始尝试将GNSS与GRACE这两种手段结合起来研究ΔTWS，从而产生比单独使用一种手段更可靠，且时空分辨率更高的ΔTWS。然而，由于两种大地测量仪器的建模方法和观测方式存在差异，如何合理解决两者空间位置和距离相关的系统偏差仍是一个悬而未决的问题。在联合反演模型中如果不能正确处理这种系统差异，意味着在反演模型中引入了不切实际的观测值，对反演结果的准确性和可靠性有极大的影响。
因此，现有技术还有待于改进和发展。
发明内容
本发明要解决的技术问题在于，针对现有技术的上述缺陷，提供一种顾及比例因子的陆地水储量变化反演方法，旨在解决现有技术中在联合反演模型中不能正确处理系统差异，导致影响了反演结果的准确性和可靠性的问题。
本发明解决技术问题所采用的技术方案如下：
第一方面，本发明提供一种顾及比例因子的陆地水储量变化反演方法，其中，所述方法包括：
获取GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据；
对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列；
根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；
根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；
根据所述GNSS观测位移和GRACE模型预测位移，得到真实GNSS观测台站位置的比例因子；
根据所述GNSS真实观测台站位置的比例因子，内插得到虚拟台站的比例因子；
根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立的联合反演ΔTWS模型；
根据所述反演ΔTWS模型，得到陆地水储量变化。
在一种实现方式中，所述对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列，包括：
将所述GNSS三维坐标时间序列的直角坐标系转换为站心坐标系，得到GNSS的NEU坐标时间序列；
采用四分位距法对所述GNSS的NEU坐标时间序列进行粗差剔除，得到剔除后的NEU坐标时间序列；
建立GNSS三维坐标时间序列模型，并根据所述GNSS三维坐标时间序列模型对所述剔除后的NEU坐标时间序列进行整理，得到处理后的GNSS三维坐标时间序列。
在一种实现方式中，所述根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移，包括：
求解变分贝叶斯独立成分分析框架中未知参数的后验概率分布，得到GNSS三维坐标时间序列的源信号概率分布函数；
根据所述源信号概率分布函数计算隐藏变量的后验概率，得到未知参数的估计值；
根据所述未知参数的估计值，得到独立成分分析模型；
通过所述独立成分分析模型模，得到水文负荷引起的GNSS观测位移。
在一种实现方式中，所述根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移，包括：
计算单位水文质量负荷在地球表面引起的径向位移和水平位移；
根据所述径向位移和水平位移得到垂直和水平负荷格林函数；
根据所述垂直和水平负荷格林函数和GRACE-Mascon产品的等效水高数据，得到所述GRACE模型预测位移。
在一种实现方式中，所述根据所述GNSS观测位移和GRACE模型预测位移，得到真实GNSS观测台站的比例因子，包括：
计算所述GNSS观测位移和GRACE模型预测位移之间的比值，将所述比值作为所述真实GNSS台站的比例因子。
在一种实现方式中，所述根据所述GNSS真实观测台站位置，内插得到虚拟台站的比例因子，包括：
根据所述GNSS真实观测台站位置生成泰森多边形；
根据所述泰森多边形和缓冲区，得到所述虚拟台站；
根据所述GNSS观测位移和GRACE模型预测位移之间的皮尔逊相关系数，得到虚拟台站的比例因子为
其中，i为邻近虚拟台站的GNSS台站个数；ρi代表第i个台站的GNSS观测和GRACE模型预测位移的PCC；βi表示第i个台站的比例因子。
在一种实现方式中，根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立的联合反演ΔTWS模型为：
其中，Av和Au为基于GNSS真实观测台站生成的垂向和水平向的格林函数矩阵；bn,be,bu分别为GNSS站心坐标系(ENU)下的所述处理后的GNSS三维坐标时间序列；bG为虚拟台站的观测值；W为GNSS观测的加权矩阵，根据GNSS观测值的中误差来定义；U是GRACE观测的加权矩阵；加权矩阵W中的元素σe，σn，σu分别为站心坐标系ENU下的不确定度；加权矩阵的U元素i表示GNSS真实观测台站的数量；βj为虚拟台站的比例因子；j表示虚拟台站的数量。
第二方面，本发明实施例还提供一种顾及比例因子的陆地水储量变化反演装置，其中，所述装置包括：
数据获取模块，用于获取GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据；
GNSS三维坐标时间序列获取模块，用于对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列；
GNSS观测位移获取模块，用于根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；
GRACE模型预测位移获取模块，用于根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；
真实GNSS观测台站比例因子获取模块，用于根据所述GNSS观测位移和GRACE模型预测位移，得到真实GNSS观测台站的比例因子；
虚拟台站比例因子获取模块，用于根据所述GNSS真实观测台站位置，内插得到虚拟台站的比例因子；
模型建立模块，用于根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立联合反演ΔTWS模型；
推演模块，用于根据所述反演ΔTWS模型，得到陆地水储量变化。
第三方面，本发明实施例还提供一种智能终端，其中，所述智能终端包括存储器、处理器及存储在所述存储器中并可在所述处理器上运行的顾及比例因子的陆地水储量变化反演程序，所述处理器执行所述顾及比例因子的陆地水储量变化反演程序时，实现如以上任一项所述的顾及比例因子的陆地水储量变化反演方法的步骤。
第四方面，本发明实施例还提供一种计算机可读存储介质，其中，所述计算机可读存储介质上存储有顾及比例因子的陆地水储量变化反演程序，所述顾及比例因子的陆地水储量变化反演程序被处理器执行时，实现如以上任一项所述的顾及比例因子的陆地水储量变化反演方法的步骤。
有益效果：与现有技术相比，本发明提供了一种顾及比例因子的陆地水储量变化反演方法，首先，获取GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据，本发明能够提供对ΔTWS的综合观测，而不是单一的水文成分。相比较传统的基于点观测的方法和基于水文模型的方法，能更全面地反映区域ΔTWS的信息，有助于区分气候变化和人类活动对水资源总量的影响。然后，根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移，以恢复由水文运动引起的地壳垂向形变，这能够保持GNSS时间序列中非水文变形源的最小污染。其次，利用负荷格林函数提取GRACE产品中的模型预测位移，并根据最小二乘确定根据所述GNSS观测位移和GRACE模型预测位移，得到两种观测手段之间的比例因子，这可以在不改变传统方法的形式的情况下引入了新的GRACE观测数据，相比较单一的数据拥有更高的空间分辨率和可靠性，可以捕捉到使用单一观测无法捕捉到的细节特征。通过校正GNSS和GRACE观测之间由于空间敏感性不同所造成的系统偏差，从而最小化两种观测手段之间的差异。最后，通过将比例因子调整过的GRACE模型预测位移引入GNSS反演模型获取ΔTWS，以生成更可靠的ΔTWS，解决了由于局部地区GNSS站点密度低而导致的不切实际的ΔTWS估计，GRACE观测的加入为稀疏GNSS台站分布的地区提供了重要的数据约束，可以在不增加新的地面观测的情况下对ΔTWS进行准确的估计。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的顾及比例因子的陆地水储量变化反演方法流程示意图。
图2是本发明实施例提供的顾及比例因子的陆地水储量变化反演方法的框架流程图。
图3是本发明实施例提供的从vbICA导出的时间函数(V)和空间响应(U)示意图。
图4是本发明实施例提供的位移和比例因子示意图。
图5是本发明实施例提供的显示了虚拟台站选址方案示意图。
图6是本发明实施例提供的显示了GNSS垂直地壳位移和GRACE模型预测位移之间的PCC和比例因子倒数的空间分布示意图。
图7是本发明实施例提供的等效水高的年振幅示意图。
图8是本发明实施例提供的显示了云南金沙江、澜沧江、南盘江流域的等效水高变化示意图。
图9是本发明实施例提供的顾及比例因子的陆地水储量变化反演装置的原理框图。
图10是本发明实施例提供的智能终端的内部结构原理框图。
具体实施方式
为使本发明的目的、技术方案及效果更加清楚、明确，以下参照附图并举实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
本技术领域技术人员可以理解，除非特意声明，这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是，本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件，但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。应该理解，当我们称元件被“连接”或“耦接”到另一元件时，它可以直接连接或耦接到其他元件，或者也可以存在中间元件。此外，这里使用的“连接”或“耦接”可以包括无线连接或无线耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的全部或任一单元和全部组合。
本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)，具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语，应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样被特定定义，否则不会用理想化或过于正式的含义来解释。
监测陆地水储量变化(ΔTWS)有助于分析极端气象水文事件发生的频率和强度变化，进而给出可能的应对和预防策略。传统基于点观测的方法(如水文观测站、测井)可以实现对不同的水成分进行独立测量，但由于其时空覆盖较差，无法在全球或区域范围内获取准确的ΔTWS。传统GNSS反演ΔTWS的方法适用于密集台站覆盖的区域，并且反演结果受到现有台站分布和密度的影响，因此在稀疏GNSS台站覆盖的场景下无法获取准确的ΔTWS。
为了解决上述技术问题，本发明提出利用比例因子校正的GRACE模型预测位移来作为额外的约束引入反演模型，以校正因地面观测不足导致的局部不切实际的反演结果，依赖于两种观测手段之间的比例因子的约束，从而可以最大程度减少观测手段之间的系统差异，提供更准确的ΔTWS，对稀疏GNSS台站覆盖区域的水文应用和水资源管理有重要的意义。
示例性方法
本实施例提供一种顾及比例因子的陆地水储量变化反演方法。
如图1所示，所述方法包括如下步骤：
步骤S100、获取GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据；
具体地，全球导航卫星系统(Global Navigation Satellite System，GNSS)，又称全球卫星导航系统，是能在地球表面或近地空间的任何地点为用户提供全天候的3维坐标和速度以及时间信息的空基无线电导航定位系统。其包括一个或多个卫星星座及其支持特定工作所需的增强系统，其利用全球卫星导航系统中载波相位差分技术(RTK)，其测量精度可以达到以厘米为单位，与传统的人工测量相比，其拥有精度高、易操作、测量设备便携、可全天候操作、测量点之间无须通视等人工测量无法比拟的优势，已广泛应用于大地测量、地壳运动、资源勘查、地籍测量等领域。
由于卫星重力观测技术的快速发展，特别是自2002年重力恢复和气候实验(GRACE)卫星的成功发射以来，地球重力场模型从根本上得到了重大改善。GRACE提供了前所未有的覆盖全球的高精度时变重力场，为地球科学中多个领域的科学研究提供了全新的观测数据和研究方法。
在本实施例中，将GNSS与GRACE这两种手段结合起来研究ΔTWS，以产生比单独使用一种手段更可靠，且时空分辨率更高的ΔTWS，需要获取GNSS和GRACE数据，具体包括GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据。
步骤S200、对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列；
在一种实现方式中，本实施例所述步骤S200包括如下步骤：
步骤S201、将所述GNSS三维坐标时间序列的直角坐标系转换为站心坐标系，得到GNSS的NEU坐标时间序列；
具体地，通过GNSS获取的是空间直角坐标系(WGS-84)下的XYZ坐标，在时间序列分析中，需要将其转化为站心坐标系ENU坐标。转换公式如式：
其中，e，n，u分别为ENU坐标系下的三个分量的坐标；x，y，z为对应的WGS-84坐标系下的坐标；x0，y0，z0分别为测站初始位置的WGS-84坐标；M为转换矩阵，具体表达式如式：
其中，λ0和φ0表示NEU坐标系下的大地经度和大地纬度。
步骤S202、采用四分位距法对所述GNSS的NEU坐标时间序列进行粗差剔除，得到剔除后的NEU坐标时间序列；
具体地，由于测站位置解算受到多路径效应、电离层延迟、轨道建模等多种因素的影响，可能会出现粗差点，本发明使用四分位距法(Interquartile Range,IQR)对粗差进行探测和剔除。IQR的基本原理如下：
IQR＝Q2-Q1
其中，将时间序列按照从小到大排列，Q1和Q2表示最靠近1/4和3/4处的下分位值和上分位值，粗差的探测区间为[Q1-1.5·IQR,Q2+1.5·IQR]，不在此区间的值被视为异常值。
步骤S203、建立GNSS三维坐标时间序列模型，并根据所述GNSS三维坐标时间序列模型对所述剔除后的NEU坐标时间序列进行整理，得到处理后的GNSS三维坐标时间序列。
具体地，对GNSS三维坐标时间序列进行建模的目的是为了剔除噪声，校正阶跃，并提取有用的信号。本实施例使用下式对时间序列进行拟合：
其中，ti表示年；a表示初始坐标；b表示线性运动速率；c和d表示周年运动振幅；e和f表示半周年运动振幅；gj表示因为更换仪器、地震等原因造成的阶跃；H(ti-Tgj)表示阶跃函数；vi为模型拟合残差。
再将剔除后的NEU坐标时间序列进行整理。
假设由N个GNSS连续站组成的区域网络，预处理后的GNSS坐标时间序列包含ENU三个分量，M＝3N为时间序列总数，T为记录的总历元数，则整理之后的数据矩阵X为：
其中，M为历元个数，T为观测台站数，x代表观测值。
举例说明，在本实施例中，对GNSS坐标时间序列进行预处理。使用德国波兹坦地学中心(GFZ)提供的CF参考框架下的非潮汐大气负荷(NTAL)与非潮汐海洋负荷(NTOL)形变产品来消除坐标时间序列中与非潮汐大气和海洋负荷相关的运动。然后，最小二乘方法被用来拟合三维坐标时间序列中的线性趋势、年度和半年周期的运动、以及地震和更换仪器产生的阶跃。最后，从原始时间序列中去除了共模误差、长期线性趋势以及阶跃，并使用剩余的GNSS时间序列来进行进一步的处理。
步骤S300、根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；
在一种实现方式中，本实施例所述步骤S300包括如下步骤：
步骤S301、求解变分贝叶斯独立成分分析框架中未知参数的后验概率分布，得到GNSS三维坐标时间序列的源信号概率分布函数；
具体地，求解变分贝叶斯独立成分分析框架中未知参数的后验概率分布。
GNSS坐标时间序列数据由几个源信号(构造形变、地球物理效应引起的非构造形变、噪声等)的线性组合得到，可以表示为：
XM*T＝AM*LSL*T+ε
其中，X为观测矩阵；A为混合矩阵；S为源信号矩阵；ε为噪声矩阵。其中变分贝叶斯独立成分分析的目的就是在源信号矩阵和混合矩阵均未知的情况下，仅通过观测矩阵X来得到源信号矩阵S和混合矩阵A。在贝叶斯网络中，贝叶斯推论用来计算未知参数的后验概率分布，而独立成分分析则是在目标函数指导下对模型隐藏变量(即源信号)进行学习。因此，二者在贝叶斯框架下是一致的。源信号是相互独立的，源矩阵S的概率分布函数可以写成：
其中，P代表概率分布函数；Si代表源信号；L表示源信号的个数。
步骤S302、根据所述源信号概率分布函数计算隐藏变量的后验概率，得到未知参数的估计值；
在独立成分分析模型中，理论上可以通过计算隐藏变量的后验概率来实现，可以用下式实现：
其中，M表示所选择的某一具体模型；p(W|M)是信号源模型，p(X|M)是规范因子，通常称为模型M的边际概率。假设模型M中所有未知参数都用向量W＝(ε,A,S,q,θ)来表示，其中，ε,A,S,q,θ分别为模型背景噪声，混合矩阵，模拟信源，混合高斯分量和模型参数，在已知观测信号X的前提下，参数W的后验概率可以表示为：
步骤S303、根据所述未知参数的估计值，得到独立成分分析模型；
步骤S304、通过所述独立成分分析模型模，得到水文负荷引起的GNSS观测位移。
Choudrey等(2002)将式(9)转为对模型负变自由能F的计算，可以下式来表示：
F[X]＝[ln(p(X,W)]p′(W)+H[X]
(10)
其中，p′(W)是后验概率p(W|X,M)的近似估计。可以进一步用因子分解的形式表示：
Miskin等(2000)证明，式(10)中的估计后验概率p′(W)最接近真实后验概率p(X|M)时，F[X]就达到极值，然后分别求p′(W)对各参数的偏导数，就可得到模型M中未知参数的估计值，通过位置参数的估计值可以求出源信号。
举例说明，在得到预处理的坐标时间序列之后，需要提取准确时间序列中的水文负荷形变。图3分别展示了5个IC的归一化时间响应、功率谱密度以及空间响应，即从vbICA导出的时间函数(V)和空间响应(U)。上方显示每个IC的时间函数和相应的功率谱密度图。时间函数在0-1之间归一化，每个IC的权重显示在左下角。下方显示空间响应，红色箭头代表水平响应，而色标代表垂直响应。空间响应代表站点运动的方向，其值与每个站点的特定分量的重要性成正比。其中IC1代表了一个共模年度信号，其空间响应显示出所有台站一致的运动趋势，主要代表了水文循环所引起的地表垂向运动。同样的，IC2和IC4也为共模信号，即台站表现出在同一方向上的周期性运动。IC2主要是NW-SE方向，IC4主要是NE-SW方向。其中，IC2的时间响应的功率谱密度表现出较大功率的长周期信号结合较小功率的年度和半年度周期信号，并且其90％以上的台站空间响应一致。所以IC2被归因于由局部小尺度水文信号，因为水平分量对于小波长异质性相比较垂直测量更加敏感。IC4的时间响应的功率谱密度在交点年误差信号周期附近出现峰值(351.5和175.9days)，因此将该分量归因于GNSS测量中的交点年误差。IC3和IC5的功率谱密度均表现出一定的年度和低频信号。IC3主要表现出较大功率的低频信号，IC5主要为年度信号。所以IC3和IC5两者之间可能存在串扰(一个IC泄漏到另外一个IC中)，但这不影响水文负载形变信号的重建。IC3和IC5的空间响应均呈现出NE-SW异向型，两者的组合记录了滇西南和滇东北的流域的小尺度水文信息。
步骤S400、根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；
在一种实现方式中，本实施例所述步骤S400包括如下步骤：
步骤S401、计算单位水文质量负荷在地球表面引起的径向位移和水平位移；
具体地，确定质量负荷在地球表面引起的径向位移。
假设将质量均匀且角半径放置在球对称、自引力、弹性且各向同性的地球表面，在地球表面引起的的径向位移由下式给出：
′
un＝Un(r)Pn(cosθ)eiωt＝ahn Pn(cosθ)eiωt
其中，un是径向位移；Un(r)是球谐展开的径向系数；Pn(cosθ)
′
为n阶勒让德多项式；hn为负荷勒夫数；eiωt表示所施加的负载的时间演变。a为地球半径；由于边界条件是根据地球质量来制定的，因此，上式又可以改写为式：
其中，mE为地球的质量；m′为定量的质量负荷，将其作为单位质量的负荷，则有：
步骤S402、根据所述径向位移和水平位移得到垂直和水平负荷格林函数；
具体地，计算垂直和水平负荷格林函数。施加在地球表面的1kg的负荷引起的垂直位移可以表示为式：
同样地，施加1kg的负荷引起地表的水平位移可以表示为式:
值得注意的是，水平位移并不适用于度数为0的情况，因此v的度数从n＝1开始。
步骤S403、根据所述垂直和水平负荷格林函数和GRACE-Mascon产品的等效水高数据，得到所述GRACE模型预测位移。
具体地，计算GRACE-Mascon产品中的质量负荷形变，根据负荷格林函数以及GRACE-Mascon产品提供的等效水高产品可以根据公式Gm＝d来计算GRACE模型预测位移，其中，G为根据研究区域划分的格网和GNSS台站生成的格林函数矩阵；m为GRACE-Mascon产品提供的研究区域的等效水高；d即为要计算的每个台站的GRACE模型预测位移。
步骤S500、根据所述GNSS观测位移和GRACE模型预测位移，得到真实GNSS观测台站位置的比例因子；
在一种实现方式中，本实施例所述步骤S500包括如下步骤：
步骤S501、计算所述GNSS观测位移和GRACE模型预测位移之间的比值，将所述比值作为所述真实GNSS观测台站位置的比例因子。
具体地，确定GRACE模型预测位移d与GNSS观测到的位移时间序列x之间的比例因子，即确定GNSS观测位移和GRACE模型预测位移之间的缩放系数，可以通过一个线型模型来实现：
x＝β·d+b
其中，β代表比例因子或缩放系数；b为平移系数。
举例说明，基于PREM模型计算负荷格林函数，并将GRACE-Mascon产品中的等效水高转化为GNSS台站位置的模型预测位移，之后将所有GNSS台站观测到的垂直地表形变整理为月平均序列，并与GRACE模型预测位移的时间序列对齐。然后，基于最小二乘原理解算得到两者的缩放因子。图4展示了15个被选中的站的经比例因子缩放前后GRACE模型预测位移和GNSS实测位移之间的对比。选定的GNSS台站观测的垂直地壳位移(黑线)、GRACE模型预测位移(绿线)和经过比例因子调整的GRACE模型预测位移(红色虚线)。每个子图下方的标签显示了GNSS和GRACE导出的位移之间的比例因子(1/β)和两者的PCC。从结果来看，基于比例因子的方法可以有效降低GRACE模型预测位移和GNSS实测数据在局部地区的差异。
步骤S600、根据所述GNSS真实观测台站位置的比例因子，内插得到虚拟台站的比例因子；
在一种实现方式中，本实施例所述步骤S600包括如下步骤：
步骤S601、根据所述GNSS真实观测台站位置生成泰森多边形；
步骤S602、根据所述泰森多边形和缓冲区，得到所述虚拟台站；
具体地，由于GRACE-Mascon产品的实际空间分辨率的限制，被选作虚拟台站的空间分辨率不应大于GRACE-Mascon产品实际空间分辨率(300km)。另外，选择的虚拟GNSS台站不能过多，以便GNSS观测能够主导反演过程。基于以上两点，本发明提出利用泰森多边形和缓冲区来确定虚拟台站的选址。泰森多边形又被称为冯洛诺伊图(Voronoi diagram)，是一组相邻点的垂直平分线组成的连续多边形，目前已经被广泛运用于台站选址。泰森多边形有3个特点：1.泰森多边形内仅有一个离散点。2.泰森多边形边上的点和其两边的离散点的距离相等。3.泰森多边形内的点到相应的离散点(台站位置)的距离最近。基于以上特性，虚拟台站的选址遵从以下步骤：
1.基于现有台站的空间分布生成泰森多边形，以确定每个台站的辐射区域；
2.为了保证反演的可靠性，尽量要求感兴趣的格网和最近的台站的距离不会超过100km。
3.未被缓冲区覆盖的泰森多边形的顶点被优先选中以作为虚拟台站的选址，并观察缓冲区是否覆盖了所有的感兴趣的格网。
步骤S603、根据所述GNSS观测位移和GRACE模型预测位移之间的皮尔逊相关系数，得到虚拟台站的比例因子为
其中，i为邻近虚拟台站的GNSS台站个数；ρi代表第i个台站的GNSS观测和GRACE模型预测位移的PCC；βi表示第i个台站的比例因子。
具体地，确定虚拟台站的比例因子。根据GEACE模型预测位移和GNSS观测位移之间的皮尔逊相关系数(PCC)，以及邻近虚拟台站的几个台站的比例因子根据下式来内插虚拟台站的比例因子：
其中，i表示邻近虚拟台站的GNSS台站个数；ρi代表第i个台站的GNSS观测和GRACE模型预测位移的PCC；βi表示第i个台站的比例因子。
举例说明，由于GRACE实际空间分辨率的限制，被选作虚拟台站的GRACE-VCD的空间分辨率不应大于GRACE的实际空间分辨率(300km)。另外，选择的虚拟GNSS台站不能过多，以便GNSS观测能够主导反演过程。基于以上两点，本发明基于泰森多边形和缓冲区来选择虚拟台站的位置。图5展示了我们提出的选址方案。图5(a)显示了基于现有GNSS台站(黑点)生成泰森多边形(黑线)，灰色框表示待估计的网格。图5(b)显示了基于现有GNSS台站生成半径为100公里的缓冲区(黄色圆圈)。图5(c)显示了添加的虚拟站(黑色矩形)生成缓冲区(绿色圆圈)。最终选择了8个虚拟台站来约束研究区的边缘区域，这样可以保证引入最少的台站而又可以保证反演结果的可靠性。
在求虚拟台站的比例因子之前，首先需要确定云南地区所有台站的GNSS实测时间序列和GRACE模型预测位移之间的相关性。PCC越高，说明GNSS受观测条件、误差等其他因素的影响程度越小，更多地代表了由水文循环引起的运动。图6显示了GNSS垂直地壳位移和GRACE模型预测位移之间的PCC和比例因子倒数的空间分布。图6(a)中的圆圈代表GNSS台站，颜色代表GNSS和GRACE位移之间的PCC。图6(b)中的圆圈代表GNSS台站，矩形代表虚拟台站，颜色代表比例因子的倒数，灰线代表基于现有GNSS站生成的泰森多边形，蓝线代表河流。图6显示两者的PCC介于0.65～0.87之间。GNSS实测时间序列和GRACE模型预测位移之间的比例因子介于1.49～2.85之间，比例因子越大说明GRACE越低估了水文活动造成形变位移。所以虚拟台站的比例因子是通过附近几个GNSS台站的实测位移和GRACE模型预测位移的相关性和比例因子共同确定的。
步骤S700、根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立联合反演ΔTWS模型；
在本实施例中，GNSS联合GRACE共同反演区域ΔTWS，即通过获取的GNSS观测位移，结合经过比例因子校正的GRACE模型预测位移，共同建立联合反演模型，包括以下步骤：
首先，建立仅GNSS反演ΔTWS模型(GNSS-only模型)。
通过GNSS网络可以准确获取地表区域三维形变特征。格林函数可以将ΔTWS映射为垂直和水平位移，以量化地表的弹性响应。因此，基于格林函数对于近场质量变化的敏感性，在台站分布合理的条件下，可以实现对局部精细负荷信号的有效捕捉。使用加权最小二乘反演算法来最小化以下失配函数：
其中，W为加权矩阵；A为格林函数矩阵，作为每个GNSS台站的三维形变与每个网格的质量变化相联系的载体；x为待估等效水高矢量；b为GNSS观测值矢量。为了避免相邻格网出现不切实际的空间变化，我们使用拉普拉斯算子来定义模型粗糙度，确保了相邻格网的质量的平滑变化。平滑因子被用来平衡模型粗糙度和数据失配之间的相对权重，合适的的值由L曲线给出。
其次，建立GNSS与GRACE联合反演ΔTWS模型(GNSS-GRACE模型)。
GNSS-only模型的敏感性会随距离增加而减弱。因此在台站覆盖不足或台站之间的空间距离较远时，反演效果会受到影响，这会导致局部质量变化不切实际。为了解决GNSS台站分布不足或者不均导致反演精度受限的问题，本发明将GRACE作为辅助观测手段，并将GRACE获得的等效水高导出的模型预测位移引入反演模型，为反演结果提供约束。新的方法原理可以表示为：
其中，W和U为加权矩阵；A和AG分别代表了基于GNSS和GRACE观测生成的格林函数矩阵；b和bG分别表示GNSS台站和虚拟台站(GRACE模型预测位移)的观测值；为平滑因子；为拉普帕斯算子。要使F(x)＝min，上式又可以进一步表示为：
其中，Av和Au为基于现有GNSS台站生成的垂向和水平向的格林函数矩阵；bn,be,bu分别为GNSS站心坐标系(ENU)下的三维坐标时间序列；bG为虚拟台站的位移；W为GNSS观测的加权矩阵，根据GNSS观测值的中误差来定义；U是GRACE观测的加权矩阵，通过GNSS与GRACE之间的缩放因子来定义。
其中，σe，σn，σu分别为站心坐标系ENU下的不确定度；i表示GNSS台站数量；βj为虚拟台站的缩放因子；j表示虚拟台站的数量。
再通过L曲线法确定反演模型中的平滑因子。反演ΔTWS是一个典型的不适定问题，因为需要估计的参数远远大于观测量。因此，本实施例使用Tikhonov正则化原理来解决这个问题，并由L曲线得出数据失配度和模型粗糙度的相对权重。
其中，为待估计的参数；表示拟合残差。则L曲线是由点组成的曲线，在此曲线上曲率最大的点可以选择近似最优的平滑因子。
为了评估本发明提出的GNSS-GRACE反演方法在实际中的表现，我们使用以GNSS-only模型的反演结果作为对比，并以GRACE-Mascon和GLDAS产品作为参考。图7显示：云南水资源分布呈从西南向东北逐渐递减的空间趋势，仅GNSS反演模型导出的EWH幅度范围约为40-300mm，GNSS-GRACE反演模型为55-320mm，GRACE为30-220mm，GLDAS为30-175mm。GNSS-only和GNSS-GRACE反演模型推断的ΔTWS的振幅明显大于GRACE和GLDAS数据，这是由于GRACE粗糙的空间分辨率以及GLDAS简化了水文模型导致的。我们注意到，由于南盘江流域的台站密度不足，GNSS-only和GNSS-GRACE反演模型推断的ΔTWS的振幅的最大偏差达到15-100mm，这意味着GNSS-only模型会大大低估ΔTWS变化。图8从左到右依次显示了云南金沙江、澜沧江、南盘江流域的等效水高变化，从上到下依次显示了来自GNSS-GRACE模型(黑线)、GNSS-Only模型(红线)、GRACE模型(绿线)和GLDAS(青线)导出的等效水高时间序列。灰色虚线代表了相应地每月气候学(多年逐月平均值)，灰色阴影表示每月等效水高赤字(每月实际值减去每月气候学)。在金沙江流域和澜沧江流域中，由GNSS-GRACE和GNSS-only模型推断的EWH变化具有较强的相关性，PCC分别为0.88和0.86。在这种GNSS台站分布合理的场景下，两种模型与GRACE和GLDAS同样具有较高的一致性。但是，南盘江流域这种GNSS台站空间分布较为稀疏的情况下，GRACE提供的空间约束变的非常重要。在这种情况下，GNSS-only和GNSS-GRACE算法导出的EWH差异较为明显，PCC仅为0.44。GNSS-only算法的EWH与GRACE和GLDAS的PCC分别为-0.13和-0.17。而在同样的条件下，GNSS-GRACE模型的反演结果与GRACE和GLDAS仍然具有较高的一致性，EWH之间的PCC分别达到了0.77和0.82。所以在稀疏台站分布的情况下引入GRACE观测十分必要，可以保证局部ΔTWS的反演不会发生不切实际的变化。
步骤S800、根据所述反演ΔTWS模型，得到陆地水储量变化。
综上，本实施例结合了GNSS和GRACE两种观测数据，相比较现有方法在空间分辨率和可靠性方面展现出优势，可以捕捉到使用单一观测手段无法捕捉到的细节特征。相比较现有方法提升了稀疏GNSS台站覆盖地区观测的可靠性和稳定性，可以在地面观测不足的情况下获取较为可靠的ΔTWS。本实施例中的方法在简单性和鲁棒性方面展现出优势，因为其在不改变传统方法的形式的情况下引入了新的观测数据。可推广性好，可迅速实现大面积的移植和应用。
示例性装置
如图9中所示，本实施例还提供一种顾及比例因子的陆地水储量变化反演装置，所述装置包括：
数据获取模块10，用于获取GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据；
GNSS三维坐标时间序列获取模块20，用于对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列；
GNSS观测位移获取模块30，用于根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；
GRACE模型预测位移获取模块40，用于根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；
真实GNSS观测台站的比例因子获取模块50，用于根据所述GNSS观测位移和GRACE模型预测位移，得到真实GNSS观测台站的比例因子；
虚拟台站的比例因子获取模块60，用于根据所述真实GNSS观测台站的比例因子，内插得到虚拟台站的比例因子；
模型建立模块70，用于根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立联合反演ΔTWS模型；
推演模块80，用于根据所述反演ΔTWS模型，得到陆地水储量变化。
在一种实现方式中，所述GNSS三维坐标时间序列获取模块20包括：
NEU坐标时间序列获取单元，用于将所述GNSS三维坐标时间序列的直角坐标系转换为站心坐标系，得到GNSS的NEU坐标时间序列；
剔除后的NEU坐标时间序列获取单元，用于采用四分位距法对所述GNSS的NEU坐标时间序列进行粗差剔除，得到剔除后的NEU坐标时间序列；
处理后的GNSS三维坐标时间序列获取单元，用于建立GNSS三维坐标时间序列模型，并根据所述GNSS三维坐标时间序列模型对所述剔除后的NEU坐标时间序列进行整理，得到处理后的GNSS三维坐标时间序列。
在一种实现方式中，本实施例所述GNSS观测位移获取模块30包括：
源信号概率分布函数获取单元，用于求解变分贝叶斯独立成分分析框架中未知参数的后验概率分布，得到GNSS三维坐标时间序列的源信号概率分布函数；
未知参数的估计值获取单元，用于根据所述源信号概率分布函数计算隐藏变量的后验概率，得到未知参数的估计值；
独立成分分析模型获取单元，用于根据所述未知参数的估计值，得到独立成分分析模型；
GNSS观测位移获取单元，用于通过所述独立成分分析模型模，得到水文负荷引起的GNSS观测位移。
在一种实现方式中，本实施例所述GRACE模型预测位移获取模块40包括：
径向位移和水平位移获取单元，用于计算单位水文质量负荷在地球表面引起的径向位移和水平位移；
垂直和水平负荷格林函数获取单元，用于根据所述径向位移和水平位移得到垂直和水平负荷格林函数；
GRACE模型预测位移获取单元，用于根据所述垂直和水平负荷格林函数和GRACE-Mascon产品的等效水高数据，得到所述GRACE模型预测位移。
在一种实现方式中，本实施例所述真实GNSS观测台站的比例因子获取模块50包括：
真实GNSS观测台站的比例因子获取单元，用于计算所述GNSS观测位移和GRACE模型预测位移之间的比值，将所述比值作为所述真实GNSS观测台站位置的比例因子。
在一种实现方式中，本实施例所述虚拟台站的比例因子获取模块60包括：
泰森多边形生成单元，用于根据所述GNSS真实观测台站位置生成泰森多边形；
虚拟台站获取单元，用于根据所述泰森多边形和缓冲区，得到所述虚拟台站；
虚拟台站的比例因子获取单元，用于根据所述GNSS观测位移和GRACE模型预测位移之间的皮尔逊相关系数，得到虚拟台站的比例因子为
其中，i为邻近虚拟台站的GNSS台站个数；ρi代表第i个台站的GNSS观测和GRACE模型预测位移的PCC；βi表示第i个台站的比例因子。
在一种实现方式中，本实施例所述模型建立模块70，用于根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立的联合反演ΔTWS模型为
其中，Av和Au为基于GNSS真实观测台站生成的垂向和水平向的格林函数矩阵；bn,be,bu分别为GNSS站心坐标系(ENU)下的所述处理后的GNSS三维坐标时间序列；bG为虚拟台站的观测值；W为GNSS观测的加权矩阵，根据GNSS观测值的中误差来定义；U是GRACE观测的加权矩阵；加权矩阵W中的元素σe，σn，σu分别为站心坐标系ENU下的不确定度；加权矩阵的U元素i表示GNSS真实观测台站的数量；βj为虚拟台站的比例因子；j表示虚拟台站的数量。
基于上述实施例，本发明还提供了一种智能终端，其原理框图可以如图10所示。该智能终端包括通过系统总线连接的处理器、存储器、网络接口、显示屏、温度传感器。其中，该智能终端的处理器用于提供计算和控制能力。该智能终端的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该智能终端的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种顾及比例因子的陆地水储量变化反演方法。该智能终端的显示屏可以是液晶显示屏或者电子墨水显示屏，该智能终端的温度传感器是预先在智能终端内部设置，用于检测内部设备的运行温度。
本领域技术人员可以理解，图10中示出的原理框图，仅仅是与本发明方案相关的部分结构的框图，并不构成对本发明方案所应用于其上的智能终端的限定，具体的智能终端以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种智能终端，智能终端包括存储器、处理器及存储在存储器中并可在处理器上运行的顾及比例因子的陆地水储量变化反演程序，处理器执行顾及比例因子的陆地水储量变化反演程序时，实现如下操作指令：
获取GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据；
对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列；
根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；
根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；
根据所述GNSS观测位移和GRACE模型预测位移，得到真实GNSS观测台站位置的比例因子；
根据所述GNSS真实观测台站位置的比例因子，内插得到虚拟台站的比例因子；
根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立联合反演ΔTWS模型；
根据所述反演ΔTWS模型，得到陆地水储量变化。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本发明所提供的各实施例中所使用的对存储器、存储、运营数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双运营数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。
综上，本发明公开了一种顾及比例因子的陆地水储量变化反演方法，所述方法包括：对GNSS三维坐标时间序列进行预处理；根据处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；根据GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；根据GNSS观测位移和GRACE模型预测位移，得到真实GNSS观测台站的比例因子；根据GNSS真实观测台站位置，得到虚拟台站的比例因子；对GNSS真实观测台站和虚拟台站进行棋盘测试，根据比例因子建立反演ΔTWS模型，进行反演。本发明依赖于虚拟台站的比例因子的约束，从而可以最大程度减少观测手段之间的系统差异，提供更准确的ΔTWS，对稀疏GNSS台站覆盖区域的水文应用和水资源管理有重要的意义。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
Claims (10)
1.一种顾及比例因子的陆地水储量变化反演方法，其特征在于，所述方法包括：
获取GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据；
对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列；
根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；
根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；
根据所述GNSS观测位移和GRACE模型预测位移，得到GNSS真实观测台站位置的比例因子；
根据所述GNSS真实观测台站位置的比例因子，内插得到虚拟台站的比例因子；
根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立联合反演ΔTWS模型；
根据所述反演ΔTWS模型，得到陆地水储量变化。
2.根据权利要求1所述的顾及比例因子的陆地水储量变化反演方法，其特征在于，所述对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列，包括：
将所述GNSS三维坐标时间序列的直角坐标系转换为站心坐标系，得到GNSS的NEU坐标时间序列；
采用四分位距法对所述GNSS的NEU坐标时间序列进行粗差剔除，得到剔除后的NEU坐标时间序列；
建立GNSS三维坐标时间序列模型，并根据所述GNSS三维坐标时间序列模型对所述剔除后的NEU坐标时间序列进行整理，得到处理后的GNSS三维坐标时间序列。
3.根据权利要求1所述的顾及比例因子的陆地水储量变化反演方法，其特征在于，所述根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移，包括：
求解变分贝叶斯独立成分分析框架中未知参数的后验概率分布，得到GNSS三维坐标时间序列的源信号概率分布函数；
根据所述源信号概率分布函数计算隐藏变量的后验概率，得到未知参数的估计值；
根据所述未知参数的估计值，得到独立成分分析模型；
通过所述独立成分分析模型模，得到水文负荷引起的GNSS观测位移。
4.根据权利要求1所述的顾及比例因子的陆地水储量变化反演方法，其特征在于，所述根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移，包括：
计算单位水文质量负荷在地球表面引起的径向位移和水平位移；
根据所述径向位移和水平位移得到垂直和水平负荷格林函数；
根据所述垂直和水平负荷格林函数和GRACE-Mascon产品的等效水高数据，得到所述GRACE模型预测位移。
5.根据权利要求1所述的顾及比例因子的陆地水储量变化反演方法，其特征在于，所述根据所述GNSS观测位移和GRACE模型预测位移，得到GNSS真实观测台站位置的比例因子，包括：
计算所述GNSS观测位移和GRACE模型预测位移之间的比值，将所述比值作为GNSS真实观测台站位置的比例因子。
6.根据权利要求1所述的顾及比例因子的陆地水储量变化反演方法，其特征在于，所述根据所述GNSS真实观测台站位置的比例因子，内插得到虚拟台站的比例因子，包括：
根据所述GNSS真实观测台站位置生成泰森多边形；
根据所述泰森多边形和缓冲区，得到所述虚拟台站；
根据所述GNSS观测位移和GRACE模型预测位移之间的皮尔逊相关系数，得到所述虚拟台站的比例因子为
其中，i为邻近虚拟台站的GNSS台站个数；ρi代表第i个台站的GNSS观测和GRACE模型预测位移的PCC；βi表示第i个台站的比例因子。
7.根据权利要求1所述的顾及比例因子的陆地水储量变化反演方法，其特征在于，根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立的反演ΔTWS模型为：
其中，Av和Au为基于GNSS真实观测台站生成的垂向和水平向的格林函数矩阵；bn,be,bu分别为GNSS站心坐标系(ENU)下的所述处理后的GNSS三维坐标时间序列；bG为虚拟台站的观测值；W为GNSS观测的加权矩阵，根据GNSS观测值的中误差来定义；U是GRACE观测的加权矩阵；加权矩阵W中的元素σe，σn，σu分别为站心坐标系ENU下的不确定度；加权矩阵的U元素i表示GNSS真实观测台站的数量；βj为虚拟台站的比例因子；j表示虚拟台站的数量。
8.一种顾及比例因子的陆地水储量变化反演装置，其特征在于，所述装置包括：
数据获取模块，用于获取GNSS三维坐标时间序列、GNSS真实观测台站位置和GRACE-Mascon产品的等效水高数据；
GNSS三维坐标时间序列获取模块，用于对所述GNSS三维坐标时间序列进行预处理，得到处理后的GNSS三维坐标时间序列；
GNSS观测位移获取模块，用于根据所述处理后的GNSS三维坐标时间序列，得到水文负荷引起的GNSS观测位移；
GRACE模型预测位移获取模块，用于根据所述GRACE-Mascon产品的等效水高数据，得到水文负荷引起的GRACE模型预测位移；
真实GNSS观测台站比例因子获取模块，用于根据所述GNSS观测位移和GRACE模型预测位移，得到真实GNSS观测台站位置的比例因子；
虚拟台站比例因子获取模块，用于根据所述GNSS真实观测台站位置的比例因子，内插得到虚拟台站的比例因子；
模型建立模块，用于根据所述GNSS真实观测台站和虚拟台站的观测值，以及所述虚拟台站的比例因子建立联合反演ΔTWS模型；
推演模块，用于根据所述反演ΔTWS模型，得到陆地水储量变化。
9.一种智能终端，其特征在于，所述智能终端包括存储器、处理器及存储在所述存储器中并可在所述处理器上运行的顾及比例因子的陆地水储量变化反演程序，所述处理器执行所述顾及比例因子的陆地水储量变化反演程序时，实现如权利要求1-7任一项所述的顾及比例因子的陆地水储量变化反演方法的步骤。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有顾及比例因子的陆地水储量变化反演程序，所述顾及比例因子的陆地水储量变化反演程序被处理器执行时，实现如权利要求1-7任一项所述的顾及比例因子的陆地水储量变化反演方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310283956.6A CN116451440A (zh) | 2023-03-22 | 2023-03-22 | 一种顾及比例因子的陆地水储量变化反演方法及智能终端 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310283956.6A CN116451440A (zh) | 2023-03-22 | 2023-03-22 | 一种顾及比例因子的陆地水储量变化反演方法及智能终端 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116451440A true CN116451440A (zh) | 2023-07-18 |
Family
ID=87129382
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310283956.6A Pending CN116451440A (zh) | 2023-03-22 | 2023-03-22 | 一种顾及比例因子的陆地水储量变化反演方法及智能终端 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116451440A (zh) |
-
2023
- 2023-03-22 CN CN202310283956.6A patent/CN116451440A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Mémin et al. | Correcting GPS measurements for non-tidal loading | |
| Brasington et al. | Monitoring and modelling morphological change in a braided gravel‐bed river using high resolution GPS‐based survey | |
| Jin et al. | Downscaling AMSR-2 soil moisture data with geographically weighted area-to-area regression kriging | |
| Bogusz et al. | Modelling the velocity field in a regular grid in the area of Poland on the basis of the velocities of European permanent stations | |
| Lovas et al. | Terrestrial laser scanning in deformation measurements of structures | |
| CN114417646B (zh) | 一种高维异构降水数据融合方法及系统 | |
| CN114004104B (zh) | 基于棋盘格测试的cors站点选取方法 | |
| Falchi et al. | Global geoid adjustment on local area for GIS applications using GNSS permanent station coordinates | |
| CN113704693B (zh) | 一种高精度的有效波高数据估计方法 | |
| CN112577470A (zh) | 一种UAV与InSAR融合监测矿区动态沉陷盆地的方法和系统 | |
| Liang et al. | Vertical surface displacement of mainland China from GPS using the multisurface function method | |
| CN116609859A (zh) | 一种气象灾害高分辨率区域模式预报系统及方法 | |
| Feng et al. | A hierarchical network densification approach for reconstruction of historical ice velocity fields in East Antarctica | |
| CN118314491A (zh) | 一种基于无人机测绘的海滩侵蚀检测和定量计算方法 | |
| CN116299466B (zh) | 一种输电通道地质形变监测方法和装置 | |
| Guillet et al. | Camera orientation, calibration and inverse perspective with uncertainties: a Bayesian method applied to area estimation from diverse photographs | |
| CN116451440A (zh) | 一种顾及比例因子的陆地水储量变化反演方法及智能终端 | |
| CN115201823B (zh) | 一种利用BDS-InSAR数据融合的地表形变监测方法 | |
| CN113916181B (zh) | 表面-内部一体化变形监测装置数据处理方法 | |
| CN113899301B (zh) | 联合gnss三维形变的区域陆地水储量变化反演方法及系统 | |
| CN113268869B (zh) | 地球表层质量变化监测方法及系统 | |
| CN111611929A (zh) | 基于LiDAR和InSAR技术的河道行洪风险点识别方法、装置、服务器及存储介质 | |
| CN117930298B (zh) | 基于卫星温度和姿态误差建模的静轨卫星定位误差修正方法及装置 | |
| CN115690341B (zh) | 全球地理高程数据90m分辨率DEM计算地形因子的校正方法 | |
| CN115047406B (zh) | 一种地空链路传播衰减区域重构方法 |
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 |