CN113868855A - 融合水位数据的地下水储量变化卫星重力正演模拟方法 - Google Patents

Abstract

本方案提供了一种融合水位数据的地下水储量变化卫星重力正演模拟方法，该方法包括：获取预设区域的边界及经纬度信息，获取预设区域的重力卫星数据和非地下水组分水储量数据；经预处理后，获得对应数据的月时间序列，并对非地下水组分水储量数据进行正演模拟信号泄露，提取月时间序列的变化趋势，通过融合水位数据的迭代正演模拟方式，获得校正泄露误差后的地下水储量变化趋势的反演值及储水系数，并基于此得到地下水储量变化月时间序列。本方案对信号进行了融合水位数据的迭代正演模拟信号校正，能够更好的减少泄露误差的影响，可以获得具有更高分辨率更高精度的区域地下水储量变化。

Description

融合水位数据的地下水储量变化卫星重力正演模拟方法

技术领域

本发明涉及重力卫星与地下水科学交叉领域，特别涉及一种融合水位数 据的地下水储量变化卫星重力正演模拟方法。

背景技术

地下水是指埋藏在地表以下各种形式的重力水，它是水资源的重要组成 部分，是一种珍贵的自然资源，也是水资源管理和水循环研究中的一个重要 参数。地下水作为目前除南极外陆地上储存最为丰富的淡水资源，在人类生 产生活中占有举足轻重的地位，与人类生活、农业灌溉、工矿企业以及市政 用水联系密切。特别是在干旱、半干旱地区，地下水往往是主要的水源，有 时在局部区域甚至是唯一可利用的水源。

然而近一个世纪以来，随着人口增长和城市化的迅速发展，对地下水的 不合理开采日益严重，全球范围内出现不同程度的地下水位下降、地面沉降、 海水入侵、岩溶塌陷等一系列地质环境问题，严重制约了社会经济发展。因 此，开展区域地下水储量变化的监测与分析工作，对揭示地下水储量变化特 征，合理利用和有效管理区域水资源以及实现地下水资源的可持续利用具有 十分重要的科学意义和社会价值。

估算地下水储量变化的传统方法主要是水井观测，根据水井观测的水位 数据与所在区域的储水系数估算地下水储量变化，但由于观测水井数量有限， 空间分布不均且受到储水系数求取方法的限制，导致这种方法对人力、物力、 财力的投入巨大且获取的地下水储量变化往往难以反映区域整体情况，存在 较大不确定性。除水井观测外，也可以基于水量平衡原理进行估算或通过水 文模型进行模拟估算地下水储量变化，但这要求有全面的水文、气象等数据 资料，估算和模拟结果的精度受到数据和模拟机制的限制。

2002年以来，GRACE(Gravity Recovery and Climate Experiment)及 GRACE-FO重力卫星(于2018年发射，其任务是接替已于2017年6月退役 的GRACE重力卫星)的出现为连续监测全球时变重力场变化信息提供了直接 观测手段，利用其监测的时变重力场信息可以反演陆地水储量变化，扣除非 地下水储量变化后即可分离出地下水储量变化。作为当前唯一能监测到地表 以下水储量变化的遥感卫星，其突破了传统地面观测在空间上的局限性，为 地下水科学提供了宝贵的数据资源。

在GRACE/GRACE-FO数据处理过程中，由于球谐展开的截断处理以及 采用滤波技术去除噪声，会引起信号泄露误差，即研究区域内的信号会泄露 到周围地区，同时周围地区的信号也可能泄露到研究区域内，造成区域之间 的真实信号相互混淆干扰，因此需对泄露误差进行校正。目前，常用的 GRACE/GRACE-FO信号校正方法主要有加法校正、尺度因子校正和迭代正 演模拟校正等方法。前两者需要利用水文模型资料正演模拟 GRACE/GRACE-FO数据的处理过程，然后将正演模拟的泄露信号作为校正 GRACE/GRACE-FO信号的依据。这两种方法严重依赖于所使用的水文模型 数据，而当前水文模型数据仍存在诸多不足，特别是在人类活动强烈的地区， 难以反映实际的水储量时空变化，用具有较大不确定性的水储量变化正演模 拟结果计算泄露误差或估计尺度因子，导致校正结果也存在较大误差。迭代正演模拟校正根据是否利用关于信号空间分布的先验信息进行校正分为约束 正演模拟与非约束正演模拟，约束正演模拟相比非约束正演模拟可以更好的 恢复研究区域内的泄露信号，但恢复信号空间分辨率的能力高度依赖于真实 信号的空间分布模式，真实信号空间分布的均匀性越小，恢复信号与真实信 号的对应性越低。而地下水储量变化的空间分布往往是不均匀的，因此迭代 正演模拟校正也很难反演出高分辨率、高精度的地下水储量变化。

发明内容

针对现有技术的不足，本发明提出一种融合水位数据的地下水储量变化 卫星重力正演模拟方法，通过联合GRACE/GRACE-FO重力卫星和水位数据 进行迭代正演模拟信号校正，减少了泄露误差的影响，从而获得具有更高分 辨率更高精度的区域地下水储量变化。具体而言，本发明提供了以下的技术 方案：

一方面，本发明提供了一种融合水位数据的地下水储量变化卫星重力正 演模拟方法，该方法包括：

步骤一、获取研究区域的边界及经纬度信息，设置背景区域的边界及经 纬度信息；基于所述背景区域的边界及经纬度信息，获取月尺度的重力卫星 球谐系数数据及非地下水组分水储量数据；

步骤二、将所述月尺度重力卫星球谐系数进行滤波处理并转化为格网数 据，得到含有泄漏误差的陆地水储量变化月时间序列；并将所述月尺度非地 下水组分水储量数据转化为阶数一致的球谐系数，进行相同的滤波处理后转 化为格网数据，得到非地下水组分水储量变化月时间序列；

步骤三、将所述陆地水储量变化月时间序列减去所述非地下水组分水储 量变化月时间序列，得到背景区域地下水储量变化月时间序列；基于所述研 究区域的边界及经纬度信息，获得水位数据，并由此得到研究区域地下水位 变化月时间序列；基于所述研究区域的边界及经纬度信息，获取背景区域地 下水储量变化中的研究区域地下水储量变化月时间序列；获取所述研究区域 外且在所述背景区域内的地下水储量变化数据，得到周围区域地下水储量变 化月时间序列；

步骤四、基于所述背景区域地下水储量变化月时间序列、所述研究区域 地下水储量变化月时间序列、所述周围区域地下水储量变化月时间序列、所 述研究区域地下水位变化月时间序列，得到背景区域地下水储量变化趋势、 研究区域地下水储量变化趋势、周围区域地下水储量变化趋势、研究区域地 下水位变化趋势；

步骤五、以任意值作为所述研究区域的储水系数模拟值，将所述储水系 数模拟值乘上对应格网的研究区域地下水位变化趋势，得到研究区域地下水 储量变化趋势模拟值，所述周围区域地下水储量变化趋势作为周围区域地下 水储量变化趋势模拟值；基于所述研究区域地下水储量变化趋势模拟值、周 围区域地下水储量变化趋势模拟值，得到研究区域正演模拟后的地下水储量 变化趋势；

计算所述研究区域正演模拟后的地下水储量变化趋势与所述研究区域地 下水储量变化趋势之间的差值；

当所述差值大于给定阈值时，根据研究区域地下水位变化趋势对应改变 所述储水系数模拟值，重新计算所述研究区域地下水储量变化趋势模拟值以 及所述差值；当所述差值小于等于给定阈值时，将所述研究区域地下水储量 变化趋势模拟值作为研究区域地下水储量变化趋势反演值，将对应的储水系 数模拟值作为研究区域地下水储量变化趋势反演值所对应的储水系数值；

步骤六、将所述研究区域地下水储量变化趋势反演值对应的储水系数， 乘上研究区域地下水位变化月时间序列，得到反演后的研究区域地下水储量 变化月时间序列。

优选的，所述非地下水组分水储量数据包括土壤含水量数据、地表水储 量数据、雪水当量数据、冰川冻土数据中的一种或其任意数量的组合。

优选的，所述步骤三中，所述研究区域地下水位变化月时间序列与所述 研究区域地下水储量变化月时间序列格网范围大小一致。

优选的，所述步骤四中，所述背景区域地下水储量变化月时间序列、所 述研究区域地下水储量变化月时间序列、所述周围区域地下水储量变化月时 间序列、所述研究区域地下水位变化月时间序列的趋势通过以下方式获得：

式中：△H(t)为月时间序列；t为时间；a为常数项；b为拟合的趋势项；A_i、T_i和

分别为振幅、周期和相位，当j＝1，为周年信号，当j＝2，为半周年信 号；ε(t)为误差。

优选的，所述步骤三中，地下水储量变化通过以下公式获得：

△GWS＝△TWS-△SMS-△SWS-△OS

式中：△TWS为陆地水储量变化，△SMS为土壤含水量变化，△SWS为 地表水储量变化，△OS为其他非地下水组分水储量变化，△GWS为地下水储 量变化；所述其他非地下水组分水储量变化包括雪水当量变化等。

本发明的方案，通过联合GRACE/GRACE-FO重力卫星和水位数据进行 迭代正演模拟信号校正，能够更好的减少泄露误差的影响，相比现有的技术 方法可以获得具有更高分辨率更高精度的区域地下水储量变化。

附图说明

图1为本发明实施例中融合水位数据的地下水储量变化卫星重力正演模 拟方法的流程图；

图2为本发明实施例中融合水位数据的迭代正演模拟算法的流程图；

图3为本发明实施例中2003-2018年间地下水储量变化和水位变化空间分 布图，其中，(a):GRACE/GRACE-FO重力卫星时变信号中地下水储量变化趋 势空间分布图,(b):黄淮海平原地下水位变化趋势空间分布图；

图4为本发明实施例中利用本发明和现有方法分别反演得到的2003-2018 年间黄淮海平原地下水储量变化空间分布图，其中：(a)为本发明反演结果； (b)为现有迭代正演模拟校正方法反演结果；

图5为本发明实施例中对本发明和现有方法反演结果的验证对比图，其 中：(a)为实测地下水位数据分布情况，(b)为用于验证的实测地下水位数据结 合综合给水度计算的2003-2018年间黄淮海平原实测地下水储量变化趋势空 间分布图，表格为本发明所述技术方法和现有迭代正演模拟校正方法与实测 地下水储量变化的对比指标结果。

具体实施方式

下面将结合本发明实施例中的图，对本发明实施例中的技术方案进行清 楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是 全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创 造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在一个具体的实施例中，本发明的技术方案可以通过以下的方法步骤来 实现。此处需要指出的是，在本实施例中，以下的各个步骤之间，并不是严 格的逻辑顺序的上、下步骤关系，仅是为了对本方案的具体实现过程以及主 要的内容进行尽量详细的阐述，例如步骤S3至S6，除了在数学计算上具有前 后计算逻辑关系的数据的对应步骤具有前后关系外，其他并行关系的数据计 算上可以是通过并行步骤的方式来实现的，并不一定以固定的上下步骤关系 来实现，以下的S1至S9的数字顺序，也不应当理解为对本方案的保护范围 的限制来理解。具体而言，结合图1所示，本方案的实现方法可以如下：

步骤S1：获取预设研究区的边界及经纬度信息，并根据预设研究区信号 可能泄露的范围设置背景区的边界及经纬度信息；根据背景区的边界及经纬 度信息获取月尺度的GRACE/GRACE-FO重力卫星数据、土壤含水量数据、 地表水储量数据等非地下水组分水储量数据；

在一个更为具体的实施方式中，例如可以将所述预设研究区的边界及经 纬度信息设置为黄淮海平原内0.5°的经纬度格网数据，所述背景区的边界及 经纬度信息取经纬度范围28°N-45°N，105°E-128°E为整个背景区。

所述重力卫星数据为GRACE科学数据系统(SDS)机构发布的CSR RL05 Level-2月重力场球谐系数，土壤含水量数据采用全球陆面数据同化系统 (GLDAS)通用陆面模型(Common Land Model，CLM)的每月土壤湿度数 据集，地表水数据采用《中国水资源公报》以及本实施例预设背景区域包含 的相关流域水资源公报发布的水库蓄水量数据，其他非地下水组分水储量数 据包括GLDAS模型模拟的雪水当量。当然，在其他情形下，其他非地下水组 分水储量数据还可以包括例如大型湖泊、河流等或能够引起陆地水储量变化 的数据，这需要根据实际数据计算需求进行调整，此处不应以本实施例列举 出的内容作为限定，本领域技术人员还可以根据具体模型计算的需求进行适 当的扩展或者选择、组合。

其中，需要说明的是：

GRACE/GRACE-FO重力卫星数据包括但不限于：CSR RL05 Level-2月 重力场球谐系数、GFZ RL05 Level-2月重力场球谐系数和JPL RL05 Level-2 月重力场球谐系数等；

土壤含水量数据可采用其他水文模型数据或实测数据；

地表水储量数据包括但不仅限于：水库蓄水量数据、大型湖泊蓄水量数 据等；

其他非地下水组分水储量数据包括但不限于：雪水当量数据、冰川冻土 数据等。

步骤S2：作为一个优选的方式，所述的GRACE/GRACE-FO重力卫星数 据使用截断至一定阶数的球谐系数，例如截断至60阶，此时需对其进行数据 处理：采用卫星激光测距数据计算的结果对C20项进行替换，并扣除一定时 间内的重力场平均数值，例如可以扣除2004年1月至2009年12月平均重力 场；对球谐系数作高斯滤波处理以降低高阶噪声的影响；将滤波后的球谐系 数转化为等效水高格网数据；在所述数据处理过程中，由于球谐系数展开阶 数的限制以及滤波处理，会引起信号泄露误差，即所述等效水高格网数据为 含有泄漏误差的陆地水储量变化月时间序列。

作为一个更为优选的实施方式，在本实施例中可采用半径为200km的高 斯滤波处理，也可采用其他半径的高斯滤波处理。当然，除此之外，在对截 断后的数据进行降噪处理时，也可以采用其他现有的滤波处理算法，此处不 再一一赘述。

S3：将所述月尺度的土壤含水量数据、地表水储量数据等非地下水组分 水储量数据转化为球谐系数并截断至与重力卫星球谐系数一致的阶数，并对 该球谐系数进行与重力卫星球谐系数同样的数据处理，此过程称为正演模拟， 得到正演模拟后的土壤含水量月时间序列、地表水储量月时间序列及非地下 水组分水储量变化月时间序列；

在本实施例中，分别对GLDAS模型的月土壤湿度数据、雪水当量数据以 及公报发布的水库蓄水量数据进行正演模拟信号泄露，分别得到了正演模拟 后的土壤含水量、地表水储量及雪水当量变化的月时间序列。

S4：将所述含有泄漏误差的陆地水储量变化月时间序列减去上一步骤得 到的正演模拟后的土壤含水量、地表水储量及其他非地下水组分水储量变化 月时间序列，从而分离出GRACE/GRACE-FO重力卫星时变信号中的地下水 储量变化月时间序列，优选的，可以直接将对应数据做差，即计算公式如下：

△GWS＝△TWS-△SMS-△SWS-△OS

式中：△TWS为陆地水储量变化，△SMS为土壤含水量变化，△SWS为 地表水储量变化，△OS为其他非地下水组分水储量变化，△GWS为地下水储 量变化；

在本实施例中，将GRACE/GRACE-FO重力卫星数据反演的含有泄漏误 差的陆地水储量变化月时间序列减去正演模拟后的土壤含水量、地表水储量 及雪水当量变化的月时间序列，得到了GRACE/GRACE-FO重力卫星时变信 号中的地下水储量变化月时间序列。

S5：基于预设研究区的边界及经纬度信息获取水位数据，并由此得到研 究区域地下水位变化月时间序列；

在本实施例中，水位数据为水井观测的实测地下水位数据，对所述实测 地下水位数据需进行插值处理并扣除所述平均重力场相对应时间段的地下水 位平均值，得到与重力卫星时变信号中地下水储量变化月时间序列格网大小 范围均一致的地下水位变化月时间序列；

上述插值处理可以采用例如克里金插值法、子区平均法、反距离加权插 值法等，也可采用现有技术中的其他插值方式；

在一个优选的实施方式中，本发明通过以下的反距离加权插值的方法进 行处理：根据近邻点的加权均值来估计未知点，以插值点与已知的样本点的 距离为权重来加权平均，离插值点越近的样本点赋予的权重越大，具体计算 公式为:

式中：

为插值点s₀处的插值结果；s_i为第i个样本点；Z(s_i)是在样本点 s_i处获得的测量值；n为插值计算时要使用的插值点周围样本点的数量；λ_i为 插值计算过程中使用的各样本点的权重，该值随着样本点与插值点之间距离 的增加而减小；d_i0是插值点s₀与各样本点s_i之间的距离；p为指数值，可通过求 均方根预测误差的最小值确定其最佳值；样本点在插值过程中所占权重的大小 受参数p值的影响，即随着样本点与插值点之间距离的增加，样本点对插值点 影响的权重按指数规律减少；在插值过程中，各样点值对插值点作用的权重 大小是成比例的，这些权重值的总和为1；

S6：基于所述预设研究区及所述背景区的边界及经纬度信息获取预设研 究区外背景区内的地下水储量变化月时间序列，也即周围区域地下水储量变 化月时间序列；

在本实施例中，考虑到黄淮海平原的地下水储量变化主要集中在研究区 内，研究区外的地下水储量变化信号泄露到研究区内的影响较小，主要考虑 研究区内地下水储量变化信号向研究区外的泄露误差，因此将周围区域地下 水储量变化设置为零。

需要说明的是，根据预设研究区的不同，周围区域地下水储量变化月时 间序列可采用水文模型数据、实测数据或公报发布数据；当周围区域地下水 储量变化数据难以获得时可直接将研究区外背景区内GRACE/GRACE-FO重 力卫星时变信号中的地下水储量变化作为周围区域地下水储量变化数据。

S7：对于所述重力卫星时变信号中地下水储量变化月时间序列、所述周 围区域地下水储量变化月时间序列和所述地下水位变化月时间序列，分解包 含在其中的线性趋势信号；在一个优选的实施方式中，可采用本发明所提供 的以下基于最小二乘拟合的谐波分析方法获得月时间序列中的变化趋势：

式中：△H(t)为月时间序列；t为时间；a为常数项；b为拟合的趋势项；A_i、 T_i和

分别为振幅、周期和相位，当j＝1，为周年信号，当j＝2，为半周年信 号；ε(t)为误差。

本实施例中分别对GRACE/GRACE-FO重力卫星时变信号中地下水储量 变化月时间序列、周围区域地下水储量变化月时间序列和地下水位变化月时 间序列进行上述的处理，也可以采用现有技术中的其他方法，计算得到 GRACE/GRACE-FO重力卫星时变信号中地下水储量变化月时间序列、周围 区域地下水储量变化月时间序列和地下水位变化月时间序列中包含的趋势信 号。如图3所示，本实施例中分别得到了2003-2018年间(a)GRACE/GRACE-FO 重力卫星时变信号中地下水储量变化趋势和(b)黄淮海平原地下水位变化趋势空间分布图。

S8：基于所述GRACE/GRACE-FO重力卫星时变信号中地下水储量变化 月时间序列、周围区域地下水储量变化月时间序列和地下水位变化月时间序 列中包含的趋势信号进行融合水位数据的迭代正演模拟；所述融合水位数据 的迭代正演模拟主要利用了储水系数、地下水位变化与地下水储量变化之间 的关系，如公式所示：

△GWS＝△GWL×Sy

式中：△GWS表示地下水储量变化，即本实施例中地下水储量变化趋势模 拟值；△GWL表示地下水位变化，即本实施例中地下水位变化月时间序列趋 势；Sy为储水系数模拟值。

结合图2所示，所述融合水位数据的迭代正演模拟具体过程如下：

(1)给定符合水文地质特性范围内的任意值作为预设研究区的每个格 网点储水系数模拟值Sy0，将所述储水系数模拟值Sy0乘上对应格网的地下水 位变化月时间序列趋势△GWL，得到预设研究区地下水储量变化趋势模拟值 Model0；将预设周围区域地下水储量变化趋势作为周围区域地下水储量变化 趋势模拟值；

(2)对(1)中预设研究区地下水储量变化趋势模拟值Model0与周围 区域地下水储量变化趋势模拟值进行正演模拟，即将格网数据进行球谐展开， 截断到60阶，并作高斯200km滤波，得到预设研究区正演模拟后的地下水 储量变化趋势Model1；

(3)计算所述预设研究区正演模拟后的地下水储量变化趋势Model1与 所述重力卫星时变信号中地下水储量变化趋势Model_GRACE对应格网之间的差 值，当所述差值大于给定阈值时，根据所述差值与所述地下水位变化月时间 序列趋势△GWL对应改变所述储水系数模拟值大小作为迭代后新的预设研究 区储水系数模拟值Sy0，然后重复上述步骤；

需要说明的是，若步骤S6中，直接将研究区外背景区内 GRACE/GRACE-FO重力卫星时变信号中的地下水储量变化作为周围区域地 下水储量变化数据，则此步骤中应同时计算周围区域正演模拟后的地下水储 量变化趋势与所述重力卫星时变信号中地下水储量变化趋势Model_GRACE对应 格网之间的差值，并将差值对应加上周围区域地下水储量变化趋势模拟值作 为迭代后新的周围区域地下水储量变化趋势模拟值；若步骤S6中周围区域地 下水储量变化数据为其他设置，则迭代过程中周围区域地下水储量变化趋势 模拟值固定不变；该重力卫星时变信号中地下水储量变化趋势Model_GRACE为 研究区外背景区内GRACE/GRACE-FO重力卫星时变信号中的地下水储量变 化。

当所述差值小于等于给定阈值时，停止迭代，将此次预设研究区地下水 储量变化趋势模拟值Model0作为预设研究区地下水储量变化趋势反演值，将 此次预设研究区储水系数模拟值Sy0作为所述预设研究区地下水储量变化趋 势反演值所对应的储水系数值。

S9：基于所述融合水位的迭代正演模拟算法，可获得所述预设研究区地 下水储量变化趋势反演值所对应的储水系数值，乘上所述地下水位变化月时 间序列即可得到反演后的预设研究区地下水储量变化月时间序列。

在本实施例中，图4(a)和图4(b)分别为利用本发明方法和现有迭代正演模 拟校正方法反演得到的2003-2018年间黄淮海平原地下水储量变化趋势空间 分布图。在预设研究区内实测地下水位数据较为密集的格网中，即在图5(a) 标记的蓝色格网中随机抽取50％的实测地下水位数据参与到本发明提出的反 演方法过程中，进而得到地下水储量变化，另外50％的实测地下水位数据结 合《华北平原地下水可持续利用调查评价》书中的综合给水度计算实测地下 水储量变化作为对本发明反演结果的独立验证。其中图5(b)为2003-2018年间 黄淮海平原实测地下水储量变化趋势空间分布图。验证结果如图5中表格所 示：在本实施例中，通过与实测地下水储量变化的对比表明，利用本发明所 述技术方法得到的反演结果其各项对比指标均优于现有迭代正演模拟校正方 法的反演结果，从而获得具有更高分辨率更高精度的地下水储量变化。

在另一个具体的实施例中，本发明的技术方案还可以通过一种融合水位 数据的地下水储量变化卫星重力正演模拟系统来实现，以配合执行上述的方 法。具体而言，该系统可以通过以下的方式来实现：

该系统包括：

数据获取模块，用于获取研究区域的边界及经纬度信息，设置背景区域 的边界及经纬度信息；基于所述背景区域的边界及经纬度信息，获取月尺度 的重力卫星数据及非地下水组分水储量数据；

数据预处理模块，用于将所述月尺度重力卫星球谐系数进行滤波处理并 转化为格网数据，得到含有泄漏误差的陆地水储量变化月时间序列；并将所 述月尺度非地下水组分水储量数据转化为阶数一致的球谐系数，进行相同的 滤波处理后转化为格网数据，得到非地下水组分水储量变化月时间序列；

时间序列数据模块，用于将所述陆地水储量变化月时间序列减去所述非 地下水组分水储量变化月时间序列，得到背景区域地下水储量变化月时间序 列；基于所述研究区域的边界及经纬度信息，获得水位数据，并由此得到研 究区域地下水位变化月时间序列；基于所述研究区域的边界及经纬度信息， 获取背景区域地下水储量变化中的研究区域地下水储量变化月时间序列；获 取所述研究区域外且在所述背景区域内的地下水储量变化数据，得到周围区 域地下水储量变化月时间序列；

趋势数据模块，用于基于所述背景区域地下水储量变化月时间序列、所 述研究区域地下水储量变化月时间序列、所述周围区域地下水储量变化月时 间序列、所述研究区域地下水位变化月时间序列，得到背景区域地下水储量 变化趋势、研究区域地下水储量变化趋势、周围区域地下水储量变化趋势、 研究区域地下水位变化趋势；

迭代正演模拟模块，用于以任意值作为所述研究区域的储水系数模拟值， 将所述储水系数模拟值乘上对应格网的研究区域地下水位变化趋势，得到研 究区域地下水储量变化趋势模拟值，所述周围区域地下水储量变化趋势作为 周围区域地下水储量变化趋势模拟值；基于所述研究区域地下水储量变化趋 势模拟值、周围区域地下水储量变化趋势模拟值，得到研究区域正演模拟后 的地下水储量变化趋势；

计算所述研究区域正演模拟后的地下水储量变化趋势与所述研究区域地 下水储量变化趋势之间的差值；

当所述差值大于给定阈值时，根据研究区域地下水位变化趋势对应改变 所述储水系数模拟值，重新计算所述研究区域地下水储量变化趋势模拟值以 及所述差值；当所述差值小于等于给定阈值时，将所述研究区域地下水储量 变化趋势模拟值作为研究区域地下水储量变化趋势反演值，将对应的储水系 数模拟值作为研究区域地下水储量变化趋势反演值所对应的储水系数值；

储量计算模块，用于将所述研究区域地下水储量变化趋势反演值对应的 储水系数，乘上研究区域地下水位变化月时间序列，得到反演后的研究区域 地下水储量变化月时间序列。

优选的，所述非地下水组分水储量数据包括土壤含水量数据、地表水储 量数据、雪水当量数据、冰川冻土数据中的一种或其任意数量的组合。

优选的，所述研究区域地下水位变化月时间序列与所述研究区域地下水 储量变化月时间序列格网范围大小一致。

优选的，所述时间序列数据模块中，地下水储量变化月时间序列通过以 下公式获得：

△GWS＝△TWS-△SMS-△SWS-△OS

式中：△TWS为陆地水储量变化，△SMS为土壤含水量变化，△SWS为 地表水储量变化，△OS为其他非地下水组分水储量变化，△GWS为地下水储 量变化；所述其他非地下水组分水储量变化包括雪水当量变化等。

优选的，

所述背景区域地下水储量变化月时间序列、所述研究区域地下水储量变 化月时间序列、所述周围区域地下水储量变化月时间序列、所述研究区域地 下水位变化月时间序列的趋势通过以下方式获得：

式中：△H(t)为月时间序列；t为时间；a为常数项；b为拟合的趋势项；A_i、 T_i和

分别为振幅、周期和相位，当j＝1，为周年信号，当j＝2，为半周年信 号；ε(t)为误差。

优选的，所述系统还包括交互模块，用于接收人工数据的输入，以及系 统中间数据、结果数据的输出。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程， 是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计 算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的 流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(RandomAccess Memory，RAM)等。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其 限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术 人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者 对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术 方案的本质脱离本发明实施例技术方案的精神和范围。

Claims (9)

1.一种融合水位数据的地下水储量变化卫星重力正演模拟方法，其特征在于，该方法包括：

步骤一、获取研究区域的边界及经纬度信息，设置背景区域的边界及经纬度信息；基于所述背景区域的边界及经纬度信息，获取月尺度的重力卫星球谐系数数据及非地下水组分水储量数据；

步骤二、将所述月尺度重力卫星球谐系数转化为格网数据，得到含有泄漏误差的陆地水储量变化月时间序列；并将所述月尺度非地下水组分水储量数据转化为阶数一致的球谐系数，转化为格网数据，得到非地下水组分水储量变化月时间序列；

步骤三、将所述陆地水储量变化月时间序列减去所述非地下水组分水储量变化月时间序列，得到背景区域地下水储量变化月时间序列；基于所述研究区域的边界及经纬度信息，获得水位数据，并由此得到研究区域地下水位变化月时间序列；基于所述研究区域的边界及经纬度信息，获取背景区域地下水储量变化中的研究区域地下水储量变化月时间序列；获取所述研究区域外且在所述背景区域内的地下水储量变化数据，得到周围区域地下水储量变化月时间序列；

步骤四、基于所述背景区域地下水储量变化月时间序列、所述研究区域地下水储量变化月时间序列、所述周围区域地下水储量变化月时间序列、所述研究区域地下水位变化月时间序列，得到背景区域地下水储量变化趋势、研究区域地下水储量变化趋势、周围区域地下水储量变化趋势、研究区域地下水位变化趋势；

步骤五、以任意值作为所述研究区域的储水系数模拟值，将所述储水系数模拟值乘上对应格网的研究区域地下水位变化趋势，得到研究区域地下水储量变化趋势模拟值，所述周围区域地下水储量变化趋势作为周围区域地下水储量变化趋势模拟值；基于所述研究区域地下水储量变化趋势模拟值、周围区域地下水储量变化趋势模拟值，得到研究区域正演模拟后的地下水储量变化趋势；

计算所述研究区域正演模拟后的地下水储量变化趋势与所述研究区域地下水储量变化趋势之间的差值；

当所述差值大于给定阈值时，根据所述差值与地下水位变化月时间序列趋势，迭代更新研究区正演模拟后的地下水储量变化趋势；当所述差值小于等于给定阈值时，将所述研究区域地下水储量变化趋势模拟值作为研究区域地下水储量变化趋势反演值，将对应的储水系数模拟值作为研究区域地下水储量变化趋势反演值所对应的储水系数值；

步骤六、将所述研究区域地下水储量变化趋势反演值对应的储水系数，乘上研究区域地下水位变化月时间序列，得到反演后的研究区域地下水储量变化月时间序列。

2.根据权利要求1所述的方法，其特征在于，所述非地下水组分水储量数据包括土壤含水量数据、地表水储量数据、雪水当量数据、冰川冻土数据中的一种或其任意数量的组合。

3.根据权利要求1所述的方法，其特征在于，所述步骤三中，所述研究区域地下水位变化月时间序列与所述研究区域地下水储量变化月时间序列格网范围大小一致。

4.根据权利要求1所述的方法，其特征在于，所述步骤四中，所述背景区域地下水储量变化月时间序列、所述研究区域地下水储量变化月时间序列、所述周围区域地下水储量变化月时间序列、所述研究区域地下水位变化月时间序列的趋势通过以下方式获得：

式中：△H(t)为月时间序列；t为时间；a为常数项；b为拟合的趋势项；A_i、T_i和

分别为振幅、周期和相位，当j＝1，为周年信号，当j＝2，为半周年信号；ε(t)为误差。

5.根据权利要求1所述的方法，其特征在于，所述步骤三中，地下水储量变化通过以下公式获得：

△GWS＝△TWS-△SMS-△SWS-△OS

式中：△TWS为陆地水储量变化，△SMS为土壤含水量变化，△SWS为地表水储量变化，△OS为其他非地下水组分水储量变化，△GWS为地下水储量变化；所述其他非地下水组分水储量变化包括雪水当量变化等。

6.根据权利要求1所述的方法，其特征在于，所述步骤五中，研究区域正演模拟后的地下水储量变化趋势通过以下方式获得：

将研究区地下水储量变化趋势模拟值、周围区域地下水储量变化趋势模拟值对应的格网数据进行球谐展开，并进行截断，再作高斯滤波，得到研究区正演模拟后的地下水储量变化趋势。

7.根据权利要求1所述的方法，其特征在于，所述步骤五中，迭代更新研究区正演模拟后的地下水储量变化趋势，通过以下方式进行：

基于储水系数模拟值、地下水位变化月时间序列趋势、地下水储量变化趋势模拟值关系，根据所述差值，调整储水系数的大小，作为迭代后新的储水系数模拟值并计算研究区地下水储量变化趋势模拟值，再将研究区地下水储量变化趋势模拟值、周围区域地下水储量变化趋势模拟值对应的格网数据进行球谐展开，并进行截断及高斯滤波，得到迭代更新后的研究区正演模拟后的地下水储量变化趋势。

8.根据权利要求7所述的方法，其特征在于，储水系数模拟值、地下水位变化月时间序列趋势、地下水储量变化趋势模拟值关系为：

△GWS＝△GWL×Sy

式中：△GWS为地下水储量变化趋势模拟值；△GWL为地下水位变化月时间序列趋势；Sy为储水系数模拟值。

9.根据权利要求1所述的方法，其特征在于，所述步骤五中，对于周围区域地下水储量变化趋势模拟值的获取：

若直接将重力卫星时变信号中的地下水储量变化作为周围区域地下水储量变化月时间序列，则同时计算正演模拟后的周围区域地下水储量变化趋势与重力卫星时变信号中地下水储量变化趋势在对应格网范围内的格网差值，将所述格网差值加上周围区域地下水储量变化趋势模拟值后，作为迭代后的新的周围区域地下水储量变化趋势模拟值；

否则，迭代过程中周围区域地下水储量变化趋势模拟值固定不变。

Images