一种反演区域地下水位变化的卫星重力协同正演模拟方法

技术领域

本发明涉及重力卫星与地下水科学交叉领域，具体而言，涉及一种反演区域地下水位变化的卫星重力协同正演模拟方法。

背景技术

相比于地表水资源的时空分布不均的特点，地下水资源具有分布广泛、水质好、出水稳定、便于获取等优点，地下水资源已成为干旱及半干旱地区的主要用水来源。但随着人类社会、经济的不断发展，人类对地下水的需求日益增长，由此出现部分区域地下水的过度开采，造成许多地区含水层亏空，形成了地下水降落漏斗。因此，准确地监测地下水的变化为日后地下水资源管理提供依据已成为值得关注的问题。

目前，对地下水的检测主要有通过地下水位监测井及水文地质参数计算地下水储量变化、利用其它水文气象数据基于水量平衡原理估算地下水储量变化，以及构建地下水数值模拟模型模拟预测地下水储量变化等几种较为传统的方法。但却存在水文模型数据的不准确，所需输入数据太多，地下水位监测井的数量、分布不均等因素，使得上述方法计算结果存在较大不确定性。尤其在地下水位监测井缺乏的地区，地下水位监测井分布不足会使得其监测结果无法排除地下水位异常波动的影响，同时也难以准确、全面的计算区域地下水储量变化情况。

自2002年美国宇航局(NASA)联合德国空间局(DLR)联合研制发射重力卫星(Gravity Recovery and Climate Experiment,GRACE)以来，重力卫星被广泛运用于反演多个地区地下水储量变化并得到了验证，突破了传统地下水监测手段在空间上的局限，展现了其在监测地下水储量变化这方面的潜力和价值。

处理GRACE/GRACE-FO数据时，球谐展开并截断及滤波这些过程会引起信号的泄露，研究区内外信号的相互泄露、干扰，最终影响反演结果。为减少信号泄露，常用的方法主要有加法校正、尺度因子校正和迭代正演模拟校正等方法。其中，加法校正和尺度因子校正主要通过先验信息正演模拟GRACE数据的处理过程，以达到校正GRACE信号的目的。这两种方法可以反演出区域地下水储量变化空间分布情况，但此类方法过于依赖水文模型，同时在人类活动强烈、干旱半干旱地区，水文模型难以反映其真实的地下水储量变化，故此类方法所反演的地下水储量变化结果与水文模型相似且存在较大误差。迭代正演模拟校正根据是否利用相关信号空间分布的先验信息可分为约束迭代正演模拟校正和非约束迭代正演模拟校正。与非约束迭代正演模拟校正相比，约束迭代正演模拟可以更好的恢复研究区内的泄露信号，但是该方法所依赖的先验信息的空间分布通常为均匀分布，而地下水储量变化的空间分布往往是不均匀的，因此约束迭代正演模拟校正也很难反演出高分辨率、高精度的地下水储量变化。这些问题造成了目前重力卫星难以获取区域高分辨率地下水位变化。

发明内容

针对现有技术的不足，本发明提供一种反演区域地下水位变化的卫星重力协同正演模拟方法，通过利用给水度信息及重力卫星数据协同反演地下水位变化，能更好的恢复信号泄露，同时能有效的增加对地下水位监测井分布少区域的地下水位变化的认知。

为达到上述目的，本发明提供了一种反演区域地下水位变化的卫星重力协同正演模拟方法，其包括：

步骤S1：获取研究区范围及其经纬度信息，根据研究区的位置判定可能受到信号泄露的影响范围，由此设置缓冲区的范围并获取其经纬度信息，其中，缓冲区的范围包括研究区的范围；

步骤S2：根据缓冲区边界及经纬度信息获取重力卫星球谐数据、除地下水组分以外的非地下水组分水储量数据、研究区给水度值以及包括研究区在内的整个缓冲区模型地下水储量数据；

步骤S3：对重力卫星球谐数据进行处理后，转换为含有泄露误差的陆地水储量变化数据；

步骤S4：对除地下水组分以外的非地下水组分水储量数据进行正演模拟处理，得到除地下水组分以外的非地下水组分水储量变化数据，并将非地下水组分水储量变化数据从步骤S3得到的含有泄露误差的陆地水储量变化数据中扣除，得到含有泄露误差的地下水储量变化数据，即为重力卫星时变信号中的地下水储量变化数据；

步骤S5：建立研究区虚拟地下水位变化值，利用虚拟地下水位变化值及研究区给水度值计算研究区地下水储量变化模拟值，并将研究区以外的缓冲区的模型地下水储量数据中水储量变化数据作为周边区地下水储量变化数据；

步骤S6：基于步骤S4得到的重力卫星时变信号中的地下水储量变化数据、步骤S5中的周边区地下水储量变化数据及研究区地下水储量变化模拟值，进行协同研究区给水度的迭代正演模拟。

在本实施例中，其中，除地下水组分以外的非地下水组分水储量数据包括土壤含水量数据、地表水储量数据、雪水当量数据、冰川冻土数据中的其中一种或至少两种的组合，其中，地表水储量数据包括：水库蓄水量数据及大型湖泊蓄水量数据。

在本实施例中，其中，步骤S3中对重力卫星球谐数据进行的处理具体为：

将重力卫星球谐数据球谐展开并进行截断，再进行滤波处理，得到处理后的球谐数据。

在本实施例中，其中，步骤S4中对除地下水组分以外的非地下水组分水储量数据进行正演模拟处理，与重力卫星球谐数据进行处理的方式相同，具体为：

先对除地下水组分以外的非地下水组分水储量数据作与重力卫星球谐数据相同的球谐展开与截断，再进行滤波处理，得到除地下水组分以外的非地下水组分水储量变化数据。

在本实施例中，其中，步骤S4中得到含有泄露误差的地下水储量变化数据的计算过程具体为：

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

式中：△GWS为含有泄露误差的地下水储量变化数据，△TWS为含有泄露误差的陆地水储量变化数据，△SMS为土壤含水量变化数据，△SWS为地表水储量变化，△OS为其他非地下水组分水储量变化，其中，所述其他非地下水组分水储量变化包括雪水当量变化与冰川质量变化。

在本实施例中，其中，步骤S5中计算地下水储量变化模拟值的具体过程为：

假设任意值作为研究区地下水位变化模拟值，将获取到的研究区给水度值乘以研究区地下水位变化模拟值，得到研究区地下水储量变化模拟值，计算公式为：

ΔGWS＝ΔGWL×μ

式中，ΔGWS研究区地下水储量变化模拟值，μ为无量纲的研究区给水度系数，ΔGWL研究区地下水位变化模拟值。

在本实施例中，其中，步骤S6具体包括：

步骤S601：将研究区地下水储量变化模拟值及周边区地下水储量变化数据进行正演处理，得到研究区正演模拟后的地下水储量变化数据；

步骤S602：计算研究区正演模拟后的地下水储量变化数据与步骤S4得到的含有泄露误差的地下水储量变化数据的差值；

步骤S603：当上述差值大于预设阈值时，根据上述差值调整研究区地下水储量变化模拟值，并重复步骤S601和S602，直至上述差值小于等于预设阈值停止迭代；

步骤S604：将此时对应的研究区地下水储量变化模拟值作为研究区当前的地下水位变化数据，并将此时对应的研究区地下水储量变化模拟值作为研究区地下水储量变化趋势反演结果。

在本实施例中，其中，步骤S601中将研究区地下水储量变化模拟值进行的正演处理具体为：

先将研究区地下水储量变化模拟值作与重力卫星球谐数据相同的球谐展开与截断，再进行滤波处理，得到研究区正演模拟后的地下水储量变化数据。

在本实施例中，其中，周边区地下水储量变化数据来源可为实测数据、水文模型数据或相关公报发布的数据；

其中，当周边区地下水储量变化数据无法从上述数据中获取到或其变化值小于一设定值时，可以将重力卫星球谐数据时变信号中对应的地下水储量变化作为周边区地下水储量变化数据，并在迭代过程中对其进行正演模拟，计算正演模拟后的周边区地下水储量变化与重力卫星球谐数据时变信号中对应的地下水储量变化之间的差值，将得到的差值加上周边区地下水储量变化模拟值后，作为新的周边区地下水储量变化数据。

在本实施例中，其中，除地下水组分以外的非地下水组分水储量数据与重力卫星球谐数据的格网范围大小一致。

本发明提供的反演区域地下水位变化的卫星重力协同正演模拟方法，与现有技术相比，通过协同重力卫星和给水度进行迭代正演模拟信号校正，能更好的恢复信号泄露，同时能有效的增加对地下水位监测井分布少区域的地下水位变化认知。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明一实施例的方法流程图；

图2为本发明一实施例中协同给水度的迭代正演模拟算法的流程示意图；

图3为华北平原2003-2016年GRACE重力卫星时变信号中的地下水储量变化趋势空间分布示意图；

图4A为利用本发明方法反演得到的华北平原2003-2016年地下水位变化趋势空间分布示意图；

图4B为利用本发明方法反演得到的华北平原2003-2016年地下水储量变化趋势空间分布示意图。

具体实施方式

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

图1为本发明一实施例的方法流程图，如图1所示，本实施例提供一种反演区域地下水位变化的卫星重力协同正演模拟方法，其包括：

步骤S1：获取研究区范围(边界)及其经纬度信息，根据研究区的位置判定可能受到信号泄露的影响范围，由此设置缓冲区的范围(边界)并获取其经纬度信息，其中，缓冲区的范围包括研究区的范围；

本实施例以华北平原作为研究区为例，设定的缓冲区的经纬度信息为108°E-128°E，30°N-45°N，此设定范围仅为举例说明，并非限定本发明实施的范围。

步骤S2：根据缓冲区边界及经纬度信息获取重力卫星球谐数据、除地下水组分以外的非地下水组分水储量数据、研究区给水度值、包括研究区在内的整个缓冲区模型地下水储量数据；

在本实施例中，其中，除地下水组分以外的非地下水组分水储量数据包括土壤含水量数据、地表水储量数据、雪水当量数据、冰川冻土数据中的其中一种或至少两种的组合，其中，地表水储量数据包括：水库蓄水量数据、大型湖泊蓄水量数据等。

优选的，选用的重力卫星球谐数据可以为GRACE科学数据系统(SDS)机构发布的CSR RL06 Level-2月重力场球谐数据，也可以为国内外其它重力卫星球谐数据，本发明不对其进行限定；土壤含水量数据来源可以为采用全球陆面数据同化系统(GLDAS)通用陆面模型(Common Land Model，CLM)土壤湿度数据集，也可采用实测数据或其它水文模型数据，本发明不对其进行限定；地表水储量数据可以选用《中国水资源公报》及缓冲区内相关流域水资源公报发布的水库蓄水量数据，但本发明并不限定与此；研究区内给水度取值例如采用《华北平原地下水可持续利用调查评价》中的综合给水度值，其选用数据可根据不同研究区采用不同的数据选取方式，本发明不对其进行限定。

优选的，模型地下水储量数据可以选用包括：WGHM模型(WaterGAP全球水文学模型，WaterGAP Global Hydrology model)中的数据、PCR-GLOBWB模型(PCRaster GlobalWater Balance，全球水文模型)中的数据等不同模型地下水储量数据，或者根据研究区的不同，选用实测地下水储量数据或水资源公报发布地下水储量变化数据等，本发明不对其进行限定。任然以华北平原为例，考虑到华北平原周边地区包括山西盆地地下水储量变化会对研究区内地下水储量变化信号产生影响，因此研究区周边缓冲区地下水储量变化数据选取WGHM模型中地下水模块中的数据。当周围区域地下水储量变化数据难以获取或变化较小时可将除研究区以外的缓冲区GRACE重力卫星时变信号中的地下水储量变化作为周边区域地下水储量变化数据。

步骤S3：对重力卫星球谐数据进行处理后，转换为含有泄露误差的陆地水储量变化数据；

在本实施例中，其中，步骤S3中对重力卫星球谐数据进行的处理具体为：

将重力卫星球谐数据球谐展开并进行截断，再进行滤波处理，得到处理后的球谐数据。

具体的，可以先将GRACE重力卫星球谐数据球谐数据展开并在一定阶数截断，例如截断至60阶，在此过程中需采用卫星激光测距数据计算的结果对C20项进行替换，并扣除一定时间段内的重力场平均数值；然后对展开后的数据作滤波处理以降低高阶噪声的影响，其中，滤波处理方法包括但不限于高斯滤波、DDK滤波等，优选的，本实施例可以以半径为300km的高斯滤波处理为例，但不限于此；最后，将滤波处理后的数据转换为等效水高，此处的结果即为含有泄露误差的陆地水储量变化。

步骤S4：对除地下水组分以外的非地下水组分水储量数据进行正演模拟处理，得到除地下水组分以外的非地下水组分水储量变化数据，并将非地下水组分水储量变化数据从步骤S3得到的含有泄露误差的陆地水储量变化数据中扣除，得到含有泄露误差的地下水储量变化数据，即为重力卫星时变信号中的地下水储量变化数据；

在本实施例中，其中，步骤S4中对除地下水组分以外的非地下水组分水储量数据进行正演模拟处理，与重力卫星球谐数据进行处理的方式相同，具体为：

先对除地下水组分以外的非地下水组分水储量数据作与重力卫星球谐数据相同的球谐展开与截断，例如也截断至60阶，再进行滤波处理，得到除地下水组分以外的非地下水组分水储量变化数据。

在本实施例中，其中，步骤S4中得到含有泄露误差的地下水储量变化数据的计算过程具体为：

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

式中：△GWS为含有泄露误差的地下水储量变化数据，△TWS为含有泄露误差的陆地水储量变化数据，△SMS为土壤含水量变化数据，△SWS为地表水储量变化，△OS为其他非地下水组分水储量变化，其中，所述其他非地下水组分水储量变化包括雪水当量变化、冰川质量变化等。

步骤S5：建立研究区虚拟地下水位变化值，利用虚拟地下水位变化值及研究区给水度值计算研究区地下水储量变化模拟值，并将研究区以外的缓冲区的模型地下水储量数据中水储量变化数据作为周边区地下水储量变化数据；

在本实施例中，其中，步骤S5中计算地下水储量变化模拟值的具体过程为：

假设任意值作为研究区地下水位变化模拟值，将获取到的研究区给水度值乘以研究区地下水位变化模拟值，得到研究区地下水储量变化模拟值，计算公式为：

ΔGWS＝ΔGWL×μ

式中，ΔGWS研究区地下水储量变化模拟值，μ为无量纲的研究区给水度系数，ΔGWL研究区地下水位变化模拟值。

步骤S6：基于步骤S4得到的重力卫星时变信号中的地下水储量变化数据、步骤S5中的周边区地下水储量变化数据及研究区地下水储量变化模拟值，进行协同研究区给水度的迭代正演模拟。

其中，协同研究区给水度的迭代正演模拟方法主要通过迭代得方法调节地下水位变化模拟数值使得由此模拟的地下水储量变化无限向所述GRACE重力卫星时变信号中地下水储量变化靠近，最终反演所述研究区内地下水位及地下水储量变化。

在本实施例中，其中，步骤S6具体包括：

步骤S601：将研究区地下水储量变化模拟值及周边区地下水储量变化数据进行正演处理，得到研究区正演模拟后的地下水储量变化数据；

在本实施例中，其中，将研究区地下水储量变化模拟值进行的正演处理具体为：

先将研究区地下水储量变化模拟值作与重力卫星球谐数据相同的球谐展开与截断，例如也截断至60阶，再进行滤波处理，得到研究区正演模拟后的地下水储量变化数据。

步骤S602：计算研究区正演模拟后的地下水储量变化数据与步骤S4得到的含有泄露误差的地下水储量变化数据的差值；

步骤S603：当上述差值大于预设阈值时，根据上述差值调整研究区地下水储量变化模拟值，并重复步骤S601和S602，直至上述差值小于等于预设阈值停止迭代；

步骤S604：将此时对应的研究区地下水储量变化模拟值作为研究区当前的地下水位变化数据，并将此时对应的研究区地下水储量变化模拟值作为研究区地下水储量变化趋势反演结果。

在本实施例中，其中，周边区地下水储量变化数据来源可为实测数据、水文模型数据或相关公报发布的数据；

其中，当周边区地下水储量变化数据无法从上述数据中获取到或其变化值小于一设定值时，可以将重力卫星球谐数据时变信号中对应的地下水储量变化作为周边区地下水储量变化数据，并在迭代过程中对其进行正演模拟，计算正演模拟后的周边区地下水储量变化与重力卫星球谐数据时变信号中对应的地下水储量变化之间的差值，将得到的差值加上周边区地下水储量变化模拟值后，作为新的周边区地下水储量变化数据。

在本实施例中，其中，除地下水组分以外的非地下水组分水储量数据与重力卫星球谐数据的格网范围大小一致。

上述实施例的方法中的迭代正演过程主要宗旨在于：通过迭代的方法不断调整地下水位模拟值从而改变研究区域地下水储量变化模拟值，以期正演后的地下水储量变化模拟值能无限接近研究区域重力卫星时变信号中的地下水储量变化，由此反演研究区域的地下水位变化及地下水储量变化。

图2为本发明一实施例中协同给水度的迭代正演模拟算法的流程示意图，如图2所示，本实施例以具体数据为例将协同给水度的迭代正演模拟算法的具体过程说明如下：

(a)设定符合现实情况的任意值作为预设研究区的每个格网的地下水位变化模拟值△GWLi，将所述地下水位变化模拟值△GWLi乘上对应格网的给水度Sy，得到预设研究区地下水储量变化模拟值GWSi；将预设研究区以外的缓冲区的模型地下水储量数据中水储量变化数据作为周边区地下水储量变化值；

(b)将(a)中预设研究区地下水储量变化模拟值GWSi与周边区地下水储量变化值进行正演模拟，即将此数据进行球谐展开，截断到60阶，并作高斯300km滤波，得到预设研究区正演模拟后的地下水储量变化GWS_after i；

(c)计算预设研究区正演模拟后的地下水储量变化GWS_after i与重力卫星时变信号中地下水储量变化GWS_GRACE作差，当所述差值大于预设阈值时，根据所述差值调整所述地下水位变化模拟值作为迭代后新的预设研究区给水度模拟值GWSi+1，之后重复上述步骤；

其中，i为迭代正演次数，i取值大于等于0。

需要说明的是，若步骤S6中除研究区以外的缓冲区的地下水储量变化数据为水文模型数据、实测数据等其它来源数据，在迭代过程中周边区地下水储量变化值保持不变；若除研究区以外的缓冲区的地下水储量变化数据难以获取，则可将GRACE重力卫星时变信号中的地下水储量变化数据作为周边区地下水储量变化值，同时在该步骤迭代正演过程中需同时计算周边区正演模拟后的地下水储量变化与GWS_GRACE之间的差值，并将差值加回周边区地下水储量变化数据中作为迭代后新的地下水储量变化模拟值。

(d)当所述差值小于给定阈值时，停止迭代，将此时的研究区地下水位变化模拟值作为当前研究区基于给水度对应的地下水位变化值，给水度乘上研究区地下水位变化模拟值即可获取研究区内反演后的地下水储量变化。

图3为华北平原2003-2016年GRACE重力卫星时变信号中的地下水储量变化趋势空间分布示意图，图4A为利用本发明方法反演得到的华北平原2003-2016年地下水位变化趋势空间分布示意图，图4B为利用本发明方法反演得到的华北平原2003-2016年地下水储量变化趋势空间分布示意图，根据图3、图4A及图4B可以看出，采用本实施例通过协同重力卫星和给水度进行迭代正演模拟信号校正，能更好的恢复信号泄露，同时能有效的增加对地下水位监测井分布少区域的地下水位变化认知。

其中，图3、图4A及图4B水储量变化趋势可通过以下公式计算：

式中，ΔH(t)为月时间序列，t为时间，α为常数项，b为拟合的趋势项；A

本发明另一实施例还提供一种反演区域地下水位变化的卫星重力协同正演模拟系统，用于实现上述实施例的方法，其包括：

数据获取及预处理模块，用于获取研究区及缓冲区的边界及经纬度信息；基于所述研究区及缓冲区域的边界及经纬度信息，获取重力卫星数据、除地下水组分以外的非地下水组分水储量数据、研究区给水度；将所述重力卫星球谐数据球谐展开、截断并进行滤波处理后转化为等效水高数据，得到含有泄漏误差的陆地水储量变化数据；并对除地下水组分以外的非地下水组分水储量数据作与重力卫星球谐数据一致的处理，即球谐展开、截断至与其相同的阶数并进行滤波处理，得到正演后非地下水组分水储量变化数据；

输入数据处理模块，用于处理将输入迭代正演模拟模块的数据；将正演后非地下水组分水储量变化从含有泄漏误差的陆地水储量变化数据中扣除，得到研究区含有泄露误差的地下水储量变化数据；基于研究区及缓冲区的边界和经纬度信息，获取除所述研究区以外缓冲区域内的地下水储量变化数据，由此获得周边区域地下水储量变化数据；

迭代正演模拟模块，用于将输入数据进行组合迭代，最终反演所述研究区域地下水储量变化；基于研究区的基本信息，假设任意值作为所述研究区域的地下水位变化模拟值，利用所述研究区给水度值和所述研究区域的地下水位变化模拟值计算研究区域地下水储量变化模拟值；基于研究区域地下水储量变化模拟值、周边区地下水储量变化数据进行正演模拟，即作球谐展开、截断并滤波的处理，得到研究区正演模拟后的地下水储量变化；计算研究区正演模拟后的地下水储量变化与研究区域含有泄露误差的地下水储量变化数据之间的差值；当所述差值大于预设阈值时，根据差值改变地下水位变化模拟值，并重新计算研究区地下水储量变化模拟值；重新对研究区地下水储量变化模拟值进行正演模拟，并计算所述差值；当所述差值小于等于预设阈值时，将所述地下水位变化模拟值作为此研究区地下水位变化反演结果，将研究区地下水位变化反演结果与给水度的乘积作为研究区地下水储量变化反演结果。

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

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

本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。

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