基于机器学习负载模型提高陆地水储量反演准确性的方法

本发明属于卫星重力学、水文学等交叉技术领域，尤其涉及一种基于机器学习负载模型提高陆地水储量反演准确性的方法。

陆地水对于工业、农业和人类生活的持续发展是一个重要的资源，然而仅占全球水资源的3.5％。在中国，陆地水储量的时空分布不均匀对人们的生存和发展带来一系列的问题。特别是在中国西南地区，水资源供需矛盾引发了一系列的自然灾害(如干旱、洪涝、水土流失等)。因此，必须监测陆地蓄水异常(TWSA)情况，以评估其潜力和长期可持续性。水团的重新分布会改变一个地区的重力场，重力恢复和气候实验(GRACE)卫星可以对重力场进行监测。

GRACE卫星是2002年美国国家航空航天局主持的地球系统科学探索项目的一部分，为大尺度TWSA反演提供了有效的测量手段。已有研究表明，GRACE卫星能以前所未有的准确度监测TWSA的趋势和季节性特征。然而，由于GRACE卫星的空间分辨率较差，较难反演出小尺度区域的TWSA。此外，由于GRACE卫星系统老化，GRACE任务自2017年终止，其后续任务GRACE-FO于2018年发射，两次任务之间存在了14个月的空窗期。因此，到一种替代方法来填补GRACE和GRACE-FO之间的空白非常重要。

地壳负载的重分布会导致地表质量发生变化，导致地壳在水平方向(N、E)和垂直方向(U)上产生复杂的形变位移。并且在垂直方向上的形变尤为突出，这种关系可以利用地壳负载模型来构建。地壳负载的形变响应可以通过多种观测方法进行测量，例如全球卫星定位系统(GPS)、干涉测量技术(InSAR)以及超长基线干涉测量(VLBI)等。

GPS观测手段具有高效、准确、全天候等优点，可以准确实时地获得站点的形变。于此同时，中国地壳运动观测网(CMONOC)已经建立了18年，可以连续观测中国区域的地壳形变，并积累了大量的观测数据。前人很多研究利用CMONOC提供的数据分析了中国典型地区的地球物理现象，例如中国南部、华北平原、四川省等。但是，GPS站点的空间分布不均极大地限制了它在分析地球物理现象中的应用。此外，在中国西南地区的GPS站点主要集中在四川省和云南省。因此，如何更好地模拟未观测区域的地壳形变成为近些年的研究热点问题。其中，未观测区域表示在1°×1°的格网内不包含GPS站点的区域。

在地表质量重分布的研究中，如何能更好地建立地表负荷(如陆地蓄水量、大气和非潮汐海洋负荷)与其形变响应之间的关系是一个常见问题。此外，GPS提供了一种独立且近实时地观测地壳形变的测量手段。鉴于GPS的众多优势，近年来许多学者利用GPS观测信息反演地球物理现象取得了令人瞩目的成就。2004年Heki利用密集GPS站点序列探究了日本的季节性质量重分布问题，该研究表明利用GPS垂向时间序列可以有效地反演地壳负载的季节性变化，并可以对GRACE起到补充作用。2012年Fu等结合GPS和GRACE对尼泊尔的垂直荷载形变进行了分析，结果表明，季节形变序列与GRACE反演结果相一致。2015年Fu等阐述了一种利用GPS序列反演地面蓄水量的新方法，并论证了这种反演结果可以用来填补GRACE和GRACE-FO之间的空白。尽管利用GPS反演地壳负载取得了丰硕的成果，但很少有研究将GPS反演负载的方法应用到GPS站点分布粗糙的地区。2021年姜中山等提出了一种新的算法，它将传统的负载模型和球面斯莱皮恩(Slepian)基函数相结合，该研究表明，稀疏GPS阵列也可以作为一种反演TWSA时空变化的方法。然而，上述研究主要关注于利用GPS序列对地壳负载进行反演，而忽略了如何加密GPS站点分布的问题。因此，在没有GPS站点分布的地区，如何更好地模拟出地壳负荷形变是一个关键问题。

本发明的技术解决问题：克服现有技术的不足，提供一种基于机器学习负载模型提高陆地水储量反演准确性的方法，旨在没有GPS站点分布的未观测区域，实现地壳负荷形变的模拟，进而实现对陆地水储量TWSA的反演，以提高TWSA反演精度。

为了解决上述技术问题，本发明公开了一种基于机器学习负载模型提高陆地水储量反演准确性的方法，包括：

将研究区按照1°×1°格网进行划分，并判断格网内是否包含GPS站点；

当确定格网内包含GPS站点时，对GPS站点的垂向时间序列进行预处理，得到真实GPS垂向形变序列；

当确定格网内不包含GPS站点时，基于随机森林法，模拟研究区所有格网内的地壳垂向形变序列，得到模拟地壳垂向形变序列；

对真实GPS垂向形变序列或模拟地壳垂向形变序列进行大气和非潮汐海洋负荷改正，得到改正后的地壳垂向形变序列；

将得到的改正后的地壳垂向形变序列作为地壳负载模型的输入数据，对陆地水储量TWSA进行反演。

在上述基于机器学习负载模型提高陆地水储量反演准确性的方法中，对GPS站点的垂向时间序列进行预处理，得到真实GPS垂向形变序列，包括：

从中国大陆构造网络CMONOC提供的中国区域GPS时间序列文件中，提取得到GPS站点对应的垂向时间序列；

对提取得到的GPS站点对应的垂向时间序列进行预处理，去处由于GPSGPS站点周围地震和接收机异常所引起的阶跃项和异常值，得到真实GPS垂向形变序列。

在上述基于机器学习负载模型提高陆地水储量反演准确性的方法中，GPS站点对应的垂向时间序列的预处理过程如下：

通过式(1)，得到GPS站点对应的垂向时间序列的一阶差分ΔUgps：

ΔUgps＝Ugps(t)-Ugps(t-1)···(1)

其中，t表示时间，Ugps(t)表示t时刻GPS站点对应的垂向时间序列，Ugps(t-1)表示t-1时刻GPS站点对应的垂向时间序列；

当ΔUgps大于5mm时，表明当前阶跃异常，通过式(2)对当前阶跃进行改正，得到真实GPS垂向形变序列Udeleap：

Udeleap＝Ugps(t)+ΔUgps···(2)

其中，t＝m,m+1,…,M，m表示出现异常阶跃的序列节点位置，M表示总的GPS站点个数。

在上述基于机器学习负载模型提高陆地水储量反演准确性的方法中，基于随机森林法，模拟研究区所有格网内的地壳垂向形变序列，得到模拟地壳垂向形变序列，包括：

将GPS站点所对应的温度和气压作为输入数据，GPS站点的垂向时间序列作为输出数据，利用随机森林法依次回归M次，得到M个格网的输出结果；

对M个格网的输出结果进行平均处理，得到研究区所有格网内的模拟地壳垂向形变序列。

在上述基于机器学习负载模型提高陆地水储量反演准确性的方法中，模拟地壳垂向形变序列表示如下：

其中，Ugrid表示模拟地壳垂向形变序列，T表示随机森林法中构造的随机的、去相关的决策树的总数，gi(x)表示随机森林法中每棵决策树的反演结果。

在上述基于机器学习负载模型提高陆地水储量反演准确性的方法中，对真实GPS垂向形变序列进行大气和非潮汐海洋负荷改正的流程下：

通过式(4)，对真实GPS垂向形变序列Udeleap进行大气和非潮汐海洋负荷改正，得到改正后的地壳垂向形变序列Uhy：

Uhy＝Udeleap-UMERRA-UOMCT···(4)

其中，UMERRA表示大气负载产生的地壳形变，UOMCT表示非海洋潮汐负载所产生的地壳形变。

在上述基于机器学习负载模型提高陆地水储量反演准确性的方法中，对模拟地壳垂向形变序列进行大气和非潮汐海洋负荷改正的流程如下：

通过式(5)，对模拟地壳垂向形变序列Ugrid进行大气和非潮汐海洋负荷改正，得到改正后的地壳垂向形变序列Uhy：

Uhy＝Ugrid-UMERRA-UOMCT···(5)

其中，UMERRA表示大气负载产生的地壳形变，UOMCT表示非海洋潮汐负载所产生的地壳形变。

在上述基于机器学习负载模型提高陆地水储量反演准确性的方法中，将得到的改正后的地壳垂向形变序列作为地壳负载模型的输入数据，对陆地水储量TWSA进行反演，包括：

根据改正后的地壳垂向形变序列Uhy，估算到每日陆地水储量变化LoadTWSA：

LoadTWSA＝((Gx-Uhy)/σ)2+β2(L(x))2→min···(6)

其中，每日陆地水储量变化LoadTWSA即为陆地水储量TWSA的反演结果，Gx表示格林函数系数矩阵，σ表示水文负载形变序列的标准差，β表示平滑因子，L(x)表示拉普拉斯算子函数。

本发明具有以下优点：

利用密集分布的全球卫星定位系统(GPS)站点可精准地反演陆地水储量异常(TWSA)。然而，当GPS测站分布少或不均匀时，较大程度地限制了基于GPS反演TWSA的应用。基于此问题，本发明公开了一种基于机器学习负载模型提高陆地水储量反演准确性的方法，可以通过增加地表垂向形变的空间覆盖以提高TWSA反演准确性。

第一，在传统GPS反演TWSA的基础上利用随机森林回归方法模拟出未知格网内的地壳负荷形变，构建了有助于提高反演TWSA准确性的新型机器学习负载反演法(MLLIM)。

第二，以中国西南地区GPS序列为数据基础，分别利用MLLIM和传统GPS反演法对中国西南地区TWSA进行反演，并与重力恢复和气候实验(GRACE)卫星和全球陆地同化系统(GLADAS)模型结果进行对比。结果表明，基于MLLIM反演结果与GRACE和GRACE-FO结果的皮尔森相关系数(PCC)分别为0.91和0.88，可决系数(R2)分别为0.76和0.65；基于MLLIM反演结果与GLDAS结果的PCC和R2分别为0.79和0.64，较传统GPS反演法在PCC和R2上平均提高8.85％和7.99％，表明基于MLLIM可提高TWSA反演的准确性。

第三，应用MLLIM反演出中国西南地区各省份的TWSA，并结合降雨数据分析了研究区内各省份的TWSA变化情况，结果表明，基于MLLIM得到的TWSA与降雨数据空间分布一致，其时空变化凸起均位于云南省西南部和广西省东南部。通过基于MLLIM反演结果与GRACE和GLDAS结果进行对比，结果表明其PCC达0.86和0.94，基于MLLIM可准确地反演出GPS测站分布稀疏地区的TWSA。

基于MLLIM准确反演的TWSA可用于填补GRACE与GRACE-FO之间的空白期。新型机器学习负载反演方法的优点为反演陆地水储量变化准确性高、计算速度快。

图1是本发明实施例中一种基于机器学习负载模型提高陆地水储量反演准确性的方法的步骤流程图；

图2是本发明实施例中一种不同质量和不同半径圆盘置于地表产生垂向形变与中心距离的关系图；

图3是本发明实施例中一种中国西南地区站点信息图；其中，3(a)为中国西南地区GPS站点分布图；3(b)为中国西南地区土地利用类型图；3(c)为未知格网(G1-G87)和GPS站点间的位置关系图；3(d)为中国西南地区在中国的位置图；

图4是本发明实施例中一种GPS站点所在格网的温度、气压和垂向形变序列示意图；其中，4(a)为SCPZ站点；4(b)为SCXJ站点；4(c)为YNMJ站点；4(d)为YNML站点；

图5是本发明实施例中一种基于MLLIM反演GPS格网地壳形变结果图；其中，5(a)为模拟序列和真实序列间的PCC；5(b)为模拟序列和真实序列间的RMSE；5(c)为模拟序列和真实序列效果图；

图6是本发明实施例中一种基于MLLIM在未知格网处地壳形变模拟结果示意图；其中，6(a)为预测格网的空间位置信息；6(b)为G36格网模拟序列效果图；6(c)为G43格网模拟序列效果图；

图7是本发明实施例中一种基于MLLIM反演中国西南地区TWSA结果示意图；其中，7(a)为0.25°圆盘TWSA结果；7(b)为周年振幅空间分布图；7(c)为中国西南地区TWSA序列图；

图8是本发明实施例中一种基于MLLIM和传统GPS反演TWSA结果图；其中，8(a)～8(c)、8(d)～8(f)分别表示MLLIM及传统GPS反演法计算得到TWSA的周年振幅、半周年振幅及序列周期；

图9是本发明实施例中一种MLLIM与GRACE结果对比效果图；其中，9(a)为拟合后的平均结果；9(b)为RACE序列的周年振幅；9(c)为GRACE序列的半周年振幅；9(d)为GRAC序列的周期；9(e)为线性散点图；

图10是本发明实施例中一种MLLIM与GLDAS结果对比效果图；其中，10(a)为拟合后的平均结果；10(b)为GLDAS序列的周年振幅；10(c)为GLDAS序列的半周年振幅；10(d)为GLDAS序列的周期；10(e)为线性散点图；

图11是本发明实施例中一种基于MLLIM得到各个省份的TWSA结果图，并结合降雨数据进行分析中国西南地区的TWSA结果图；其中，11(a)～11(e)为各省份降雨空间分布图；11(f)～11(j)为结合降雨分析TWSA变化情况图。

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明公开的实施方式作进一步详细描述。

本发明的核心思想之一在于：基于机器学习和地壳负载模型，构建一种新型机器学习负载反演法(MLLIM)，并反演得到中国西南地区2011-2019年的TWSA。MLLIM主要流程如下：第一，判断1°×1°格网内是否存在GPS站点；第二，对不包含GPS站点的格网利用随机森林(RF)反演了未观测网格内的地壳形变序列；第三，利用所有地壳序列反演得到中国西南地区TWSA；第四，与GRACE和水文模型结果进行对比，验证了该MLLIM方法的准确性；第五，应用MLLIMMLLIM反演出中国西南地区各个省份的TWSA，并结合降雨数据进行分析各省份的TWSA变化关系，对西南地区水资源的有效管理和人们生活具有重要意义。

2014年，Argus等提出了利用GPS反演水储量变化的思路，不同于井位深度测量、水文模型、卫星测高及GRACE重力卫星等水储量变化监测手段。利用GPS反演陆地水储量变化应用了地壳负荷形变相关理论，可以大大增加GPS数据的物理意义，且较传统GRACE重力卫星法可有效地提高时空分辨率。然而，由于GPS站点分布不均匀，导致在GPS站点稀疏的地区反演TWSA的准确性低，GPS站点的疏密程度极大地限制了GPS反演TWSA的应用。因此，如何精准地模拟出GPS站点稀疏地区的垂向位移形变序列是提高反演精度的关键。

如图1，在本实施例中，该基于机器学习负载模型提高陆地水储量反演准确性的方法，包括：

步骤101，将研究区按照1°×1°格网进行划分。

步骤102，判断格网内是否包含GPS站点。

在本实施例中，当确定格网内包含GPS站点时，跳转执行步骤103；当确定格网内不包含GPS站点时，跳转执行步骤104。

步骤103，对GPS站点的垂向时间序列进行预处理，得到真实GPS垂向形变序列。

在本实施例中，可以从中国大陆构造网络CMONOC提供的中国区域GPS时间序列文件中，提取得到GPS站点对应的垂向时间序列。然而，由于GPS站点周围地震和接收机异常的影响，在时间序列中会存在阶跃项和异常值。因此，可以对提取得到的GPS站点对应的垂向时间序列进行预处理，去处由于GPSGPS站点周围地震和接收机异常所引起的阶跃项和异常值，进而得到真实GPS垂向形变序列。

优选的，GPS站点对应的垂向时间序列的预处理过程可以如下：

通过式(1)，得到GPS站点对应的垂向时间序列的一阶差分ΔUgps：

ΔUgps＝Ugps(t)-Ugps(t-1)···(1)

其中，t表示时间，Ugps(t)表示t时刻GPS站点对应的垂向时间序列，Ugps(t-1)表示t-1时刻GPS站点对应的垂向时间序列。

当ΔUgps大于5mm时，表明当前阶跃异常，通过式(2)对当前阶跃进行改正，得到真实GPS垂向形变序列Udeleap：

Udeleap＝Ugps(t)+ΔUgps···(2)

其中，t＝m,m+1,…,M，m表示出现异常阶跃的序列节点位置，M表示总的GPS站点个数。

步骤104，基于随机森林法(RF)，模拟研究区所有格网内的地壳垂向形变序列，得到模拟地壳垂向形变序列。

随机森林法(RF)是2001年由Breiman提出的一种用于回归、分类等任务的学习方法。首先，构造大量随机的、去相关的决策树；其次，求决策树的平均值。对于回归任务，主要优点包括：(1)基于数据可用性和用户需求对预测因子可进行选择；(2)可能包含连续分类的预测因子；(3)必须由用户指定相对较少的模型参数；(4)过度拟合的最小风险；(5)自动计算可变重要性得分，以评估单个预测因素对最终模型的贡献。

在本实施例中，可以将GPS站点所对应的温度和气压作为输入数据，GPS站点的垂向时间序列作为输出数据，利用随机森林法依次回归M次，得到M个格网的输出结果；然后，对M个格网的输出结果进行平均处理，得到研究区所有格网内的模拟地壳垂向形变序列。

优选的，模拟地壳垂向形变序列表示如下：

其中，Ugrid表示模拟地壳垂向形变序列，T表示随机森林法中构造的随机的、去相关的决策树的总数，gi(x)表示随机森林法中每棵决策树的反演结果。

步骤105，对步骤103得到的真实GPS垂向形变序列或步骤104得到的模拟地壳垂向形变序列，进行大气和非潮汐海洋负荷改正，得到改正后的地壳垂向形变序列。

在本实施例中，考虑到地壳形变序列(真实GPS垂向形变序列、模拟地壳垂向形变序列)中，有部分成分是由大气和非潮汐海洋负载所影响的，因此，为了准确地提取出水文负荷对地壳形变的影响，可以利用MERRA和OMCT模型对所有格网内的地壳形变序列进行大气和非潮汐海洋负荷改正。

优选的，对真实GPS垂向形变序列有：

通过式(4)，对真实GPS垂向形变序列Udeleap进行大气和非潮汐海洋负荷改正，得到改正后的地壳垂向形变序列Uhy：

Uhy＝Udeleap-UMERRA-UOMCT···(4)

优选的，对模拟地壳垂向形变序列有：

通过式(5)，对模拟地壳垂向形变序列Ugrid进行大气和非潮汐海洋负荷改正，得到改正后的地壳垂向形变序列Uhy：

Uhy＝Ugrid-UMERRA-UOMCT···(5)

其中，UMERRA表示大气负载产生的地壳形变，UOMCT表示非海洋潮汐负载所产生的地壳形变。

步骤106，将得到的改正后的地壳垂向形变序列作为地壳负载模型的输入数据，对陆地水储量TWSA进行反演。

地球是一个有弹性的球体，当地球表面的负载(例如：地表水、雪、冰等)发生变化时，地壳会产生相应程度的形变，这种形变是负荷形变。格林函数可用来建立负荷质量与形变之间的关系。地壳负荷形变主要表现在水平和垂直方向，然而地壳负荷形变在垂直方向更为敏感，其负荷形变振幅约为水平方向的2～3倍。格林函数Ugreen对地壳垂向负荷形变的描述如下：

其中，θ表示距离圆盘中心的角半径，Pn表示勒让德多项式，G表示牛顿万有引力常数，R表示地球半径，hn表示负荷勒夫数，g表示重力加速度。

Γn函数的推导如下：

当n＝0时，Γn函数的的表达式如下：

如图2，表示不同质量、不同半径和厚度的圆盘置于地表产生负荷形变与距离的关系图。据图2可知，近场中负荷响应非常明显。当距圆盘中心距离等于圆盘半径时，此时负荷响应仅为圆盘中心响应的1/2；当距圆盘中心距离为圆盘半径5倍时，这种响应可忽略不计。由此可见，GPS站点仅能反演出有限范围内的陆地水储量变化。因此，对于模拟未知格网的垂向负荷形变序列非常重要。

MLLIM方法是以所有改正后的地壳垂向形变序列为数据基础估算对应0.25°×0.25°格网的水储量变化。接着，利用曲率平滑算法对所获得的解进行正则化处理，并作为一组约束条件附加在求解矩阵中。具体地讲，在这个研究的每段时间内，根据改正后的地壳垂向形变序列Uhy，对抑制最小二乘问题进行最小化处理，估算到每日陆地水储量变化LoadTWSA：

LoadTWSA＝((Gx-Uhy)/σ)2+β2(L(x))2→min···(6)

其中，每日陆地水储量变化LoadTWSA即为陆地水储量TWSA的反演结果，Gx表示格林函数系数矩阵，σ表示水文负载形变序列的标准差，β表示平滑因子，L(x)表示拉普拉斯算子函数。

在上述实施例的基础上，下面对该基于机器学习负载模型提高陆地水储量反演准确性的方法(MLLIM方法)的验证进行说明。

数据和模型

本发明利用中国大陆环境监测网络提供的GPS站点的坐标时间序列进行分析。中国西南地区共有57个站点，其空间分布如图2所示。该坐标时间序列的获得是以原始观测文件、导航文件及精密星历文件为数据基础，利用GAMIT软件对其进行电离层改正、绝对天线相位中心校正、及海洋潮汐改正等，解算出站点间的基线，并用GLOBK平差软件解算出ITRF2008框架下的站点坐标，进而得到GPS站点坐标时间序列。接着对获得的垂向坐标时间序列进行阶跃项改正并剔除大于三倍标准差的序列异常值。为了避免数据缺失对TWSA反演造成的影响，本发明利用GPS序列间共有最长的时间跨度为2011-2019年。据图3可知，中国西南地区中云南省GPS站点最为密集，而在四川东部、重庆、广西等地GPS站点比较稀疏。

重力卫星数据及陆地水储量变化反演模型

利用GRACE和GRACE-FO时变重力场可以有效地反演全球区域内的陆地水储量变化，计算公式如下：

其中，Δh表示陆地水储量变化等效水高；A表示地球半径6371.39km；ρe表示地球平均密度5.51×103kg/m3；ρw表示水的密度103kg/m3；hl和kl表示l阶负荷勒夫数；Wl表示高斯平滑核函数；表示完全规格化缔合勒让德函数；ΔClm与ΔSlm表示地球重力场的球谐系数变化量。

本发明利用美国德克萨斯大学空间中心(CSR)和美国宇航局喷气推进实验室(JPL)提供的2011-2019年间的GRACE Mascon(GRACE-M)数据产品，并根据边界文件提取出中国西南地区的TWSA。为了削弱解算策略对数据产品不确定性的影响，本发明取两种数据的平均值作为最终的TWSA结果：

GLDAS数据和气压数据

为了有效地回归出未知区域的GPS垂向时间序列，本发明以全球陆地同化系统(GLADAS)提供的V2.2全球陆面同化数据和欧洲中期天气预报中心(ECMWF)提供的日尺度气压数据作为输入数据，接着分别提取出中国西南地区的日尺度地表温度和气压。其中，GLDAS V2.2陆面同化数据集的时间分辨率为天尺度，空间分辨率为0.25°，时间跨度为2003年到现在，该数据集包括温度在内的28个变量；ECMWF提供的气压数据的时间分辨率为天尺度，空间分辨率为0.1°，时间跨度为2000年到现在。图3描述了GPS站点垂向时间序列、温度时间序列及气压时间序列，以SCPZ、SCXJ、YNMJ和YNML为例。

据图4可知，温度、气压和垂向形变序列的周期长度均接近于一年。温度序列与GPS垂向形变序列成负相关的关系，且可以明显地看出气压序列与地壳垂向序列存在一定的相位差。整体而言，温度和气压时间序列的周期性较地壳形变序列更为明显，是由于影响地壳形变的因素更为复杂。

评定指标

本发明利用均方根误差(RMSE)、皮尔森相关系数(PCC)和R平方(R2)对反演结果进行评估：

其中，Y和S分别表示真实数据和模拟数据，和表示数据的平均值。RMSE描述了实际序列和模拟序列之间的离散度，如果RMSE值较小，则表明模拟结果在振幅上是稳定的。PCC是一个线性相关系数，它反映了两个量之间的线性相关程度。PCC值在-1和1之间，绝对值越接近1，相关性越强，反应其准确性越高。它主要提取序列间的季节关系。同时，R2也称为决定系数，它表示总离差平方和中可以由回归平方和解释的比值。R2值介于0和1之间，其值越大，拟合效果越好，表明准确性越高。

基于MLLIM计算TWSA

地壳形变主要可以划分为水平(N、E)方向和垂直方向(U)，其中在水平方向的序列主要体现为线性变化的特征，而地壳垂向的特征主要表现在非线性的季节性变化。GPS接收机固定于基岩上，由于热胀冷缩现象，当地表温度变化会使地表附近的基岩产生周期性的上下震动。并且由于站点海拔的不同，所在区域的气压也不同，气压的不同会影响GPS站点振幅的大小。因此，本发明利用了GLDAS V2.2提供的地表温度数据和CFWCM提供的气压数据作为输入序列，将GPS站点垂向位移序列作为输出数据进行回归。为了削弱单一站点回归误差的影响，本发明利用所有站点多次回归后平均的策略，对未知格网的垂向形变位移进行回归。

为了验证该回归方法的可用性，57个GPS站点所在格网的温度序列、气压序列和时间被作为输入数据，GPS垂向序列作为输入数据，对每个GPS垂向时间序列进行回归56次(除去该站点本身)。如图5(a)和图5(b)所示，利用PCC和RMSE两种指标对回归结果进行评估；为了整体表现该回归方法的有效性，求出了57个站点垂向模拟序列和真实序列的均值，如图5(c)所示。

据图5可知，利用该方法可以有效地对地壳负荷形变序列进行模拟。对图5(a)中的数据进一步统计，其中相关系数达0.5以上的序列占总体的84.21％，最高值高达0.79。由图5(b)中可以看出，模拟序列与真实序列间的RMSE表现良好，其中68.42％站点的RMSE低于5mm，其RMSE的值与GPS序列的质量相关。如图5(c)所示，为了综合表现该策略模拟地壳负荷形变的适用性，分别求出研究区内57个模拟和真实值序列的均值，可以看出该方法可以有效地模拟出序列周期项与振幅，其相关性达0.75，且RMSE为3.45mm。该精度可以有效地模拟出GPS垂向形变序列，因此利用该方法回归出了未知区域内的地壳垂向形变序列。每个格网回归了57次，并求其均值作为格网最终的垂向形变序列，为下一步陆地水储量反演提供了坚实的数据基础。

为了利用GPS垂向形变序列反演中国西南地区TWSA，本发明利用了西南地区改正后的57个GPS站点垂向序列和87个未知区域的模拟地壳垂向形变序列反演该地区的TWSA，其中模拟的站点点按照未知区域格网的中心点进行限制，如图6所示。

据图6可知，基于该策略可有效地模拟出未知格网垂向形变序列的周年项和半周年项。通过上述回归方法已计算出研究区内的所有格网的地壳垂向形变序列，为TWSA的反演做好准备。首先，计算出研究区内点负荷的格林函数，选取0.25°作为格网分辨率，其所对应圆盘的半径为13.90km，选取拓展边界范围2°及β＝0.01的预设参数。利用公式(5)对中国西南地区2011～2019年的陆地水储量进行反演，得到0.25°空间分辨率的陆地水储量变化值，并计算得到中国西南地区的TWSA及其变化的周年振幅，如图7所示。

据图7可知，利用新型机器学习负载反演法可反演出研究区内及周围的陆地水储量变化，图7(b)表明云南省范围内的陆地水储量周年振幅明显高于周围省份，其振幅达120mm，图7(c)表明基于新型机器学习负载反演法可反演出陆地水储量异常序列的周年项及振幅。并且云南省范围内的陆地水储量值明显大于周围的省份，其陆地水储量等效水高达120mm。为了证明反演结果的准确性，将通过该方法反演得到的陆地水储量变化分别与传统GPS反演TWSA方法、GRACE反演陆地水储量变化结果及水文模型进行对比。

MLLIM反演结果与传统GPS反演方法结果对比

为了探讨新型机器学习负载反演法与传统GPS反演方法的区别，本发明以扣除大气负荷形变的地壳垂向形变序列作为输入数据，并结合格林负荷函数，求解得到2011～2019年间中国西南地区(0.25°×0.25°)的TWSA。接着求解出各个格点序列的周年振幅、半周年振幅及序列周期，并与本发明实验结果进行对比，如图8所示。

据图8可知，基于新型机器学习负载反演法的结果与传统GPS反演TWSA结果存在明显差异。可以将图8中的子图分为三个对比组：(a)/(d)、(b)/(e)及(c)/(f)，分别表示周年振幅(I)、半周年振幅(II)及序列周期(III)三个序列特征参数。从图8组(I)部分可知，由于中国广西省境内的GPS站点较其他省份少，利用传统GPS方法反演TWSA时，由于平滑处理造成忽略一些特征点，然而在图(a)中可以明显体现出其周年振幅特征。并且这种现象在图8组(II)中体现的更加明显，可以有效地抑制由于插值造成的信号特征缺失。并分析序列的周期时(图8组III)，对于广西省和贵州省的部分存在明显的差异，利用新型机器学习负载反演法可以有效反演得到广西省周年项的漏斗情况。为了进一步判断两种反演方法的准确性，本发明将两种方法求得的TWSA结果与GRACE重力卫星反演结果和水文模型进行对比。

MLLIM反演结果与GRACE结果对比

GRACE和GRACE-FO卫星通过监测地球时变重力场的变化，能够较高精度地反演陆地水储量变化。为了验证反演结果的准确性，本发明利用JPL和CSR机构发布的Mascon数据结果，并求其均值作为最后西南地区TWSA结果。由于GRACE卫星与GRACE-FO卫星间存在空档期，因此将GRACE-M产品分成GRACE和GRACE-FO两个时段。并且GRACE反演结果的时间分辨率为月尺度，为方便与GRACE反演结果进行对比，将MLLIM反演结果进行月平均处理。

据图9可知，GRACE-M数据的计算结果在周年振幅和半周年振幅上与MLLIM反演结果相一致。其中云南省整体周年振幅最大，最高可达120mm左右，且广西省的周年振幅也较为突出，由此可见新型机器学习负载反演法对于广西省周年振幅信号(图8(a))的模拟是准确的，然而在传统GPS反演方法中就忽略了此信号。由图9(b)可以看出对于序列半周年信号的探测效果更为突出，可明显看到位于云南省、贵州省和四川中部具有半周年振幅信号的凸起。然而，由于GPS反演结果与GRACE-M数据的分辨率不同，导致在计算序列周年相位上有些差异，但整体量级均保持一致。分别求出中国西南地区两种方法的TWSA，并对序列进行拟合处理，求得GPS反演结果与GRACE和GRACE-FO数据的相关系数分别为0.91和0.88，求其R2值分别为0.76和0.65。而通过传统GPS反演法得到TWSA的结果与GRACE和GRACE-FO的PCC分别为0.87和0.79；R2值分别为0.71和0.58。因此，较传统GPS反演TWSA方法，基于MLLIM在PCC和R2的指标上平均提高7.98％和9.30％。

MLLIM反演结果与GLDAS结果对比

为了更好地验证新型机器学习负载反演法的准确性，本发明利用GLDAS V2.2数据集中的TWSA变量进行对比。其中GLDAS V2.2数据集的空间分辨率为0.5°×0.5°，其时间分辨率为天尺度，可以更好地与MLLIM反演结果进行对比。

据图10可知，GLDAS数据在周年振幅、半周年振幅和周年相位上与MLLIM反演结果整体保持一致。据图10(a)可以看出，在云南省、四川西部和广西东部的周年振幅均有明显凸起，与MLLIM反演结果(图8(a))保持一致。图10(b)中可以看出与本结果(图8(b))更为一致，均探测出了云南省、重庆市及贵州省的半周年振幅凸起点。据图10(c)可以看出位于广西省和四川东部地区周年振幅整体偏高，与本实验结果相符，进一步证明了本反演结果的周年相位的准确性。并计算出MLLIM反演结果与GLDAS数据间的PCC及R2分别为0.79和0.64。而传统GPS反演方法与GLDAS数据间的PCC及R2分别为0.72和0.60，所以基于MLLIM较传统GPS反演法在PCC及R2分别提高了9.72％和6.67％，表明MLLIM较传统GPS反演法可提高反演TWSA的准确性。

应用

本发明反演结果与GRACE和GLDAS的PCC分别为0.88和0.79，并能够有效地探测出TWSA时间序列的周年振幅、半周年振幅和周年相位信号。由于地壳垂向的负荷形变主要与陆地水储量有关，当水量增多时，地壳产生向下的垂向形变；反之，地壳产生向上的回弹形变。为了分析中国西南地区水储量变化的原因，本发明利用中国气象数据网(CMA)提供的0.5°网格产品，分别提取出每个省的平均降雨量，并与本发明反演结果、GRACE Mascon和GLDAS数据进行对比分析，如图11所示。

据图11可知，利用MLLIM反演结果与GRACE、GLDAS和降雨数据具有较好的时空一致性。其中，由于重庆市和贵州省的GLDAS数据集存在异常情况(图10(c)中可以看出)，所以将这两个地区的GLDAS数据未作对比。据图11(a)-(e)可以看出：在云南省西南部、四川中部和广西东南部降雨存在明显的振幅凸起点，这与本发明反演TWSA结果(图8(a))相一致。据图11(f)–(j)中可以看出：(1)本发明结果与GRACE、GLDAS的周期性几乎保持一致，且本发明实验结果与GLDAS数据相位保持较好，然而由于观测手段的不同，与GRACE数据存在2个月左右的相位差；(2)GPS反演结果的振幅略大于GRACE和GLDAS，主要是因为已扣除大气负荷和非潮汐海洋负荷的信号中仍具有由于观测手段所引起的共模误差等影响；(3)通过进一步的统计，得到反演结果与GRACE的相关性最大值为0.86(云南省)、最小值为0.63(四川省)，均值为0.77，与GLDAS的相关性最高为0.94(云南省)、最小值为0.78(广西省)，均值为0.86。因此，基于MLLIM可有效地反演出中国西南地区的TWSA变化情况。

然而，由于GPS观测手段的特殊性，导致GPS观测的地壳垂向形变序列中存在各种各样的噪声成分，使反演得到的TWSA存在噪声部分，因此在后续的研究会更好地去除GPS序列的噪声部分；并且由于GPS反演TWSA方式的特性，与GRACE存在一定的相位差，因此在后续的研究中会通过信号处理相位差部分，以达到与GRACE更好地拟合。

本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。