联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法

本发明涉及一种含水层组储水系数反演的方法，具体而言，涉及一种联合mt-insar技术和mgwr模型的分层含水层组储水系数反演方法，特别是针对主要因为地下水开采导致的地面沉降区，联合mt-insar技术和空间统计分析方法，量化分层含水层水位变化与地面沉降相关性，根据分层含水层水位变化对地面沉降的贡献度，将沉降数据进行分层，并结合地下水位监测数据，在空间尺度反演不同层位含水层组的储水系数。

背景技术：

1、地面沉降是指在一定的地表面积内所发生的地表海拔标高降低的现象，是一种缓变性的地质灾害，主要是由于过量开采地下水，导致地下水位降低，从而产生地面沉降的现象，造成含水层储水能力损失。不均匀的地面沉降会影响和制约当地经济建设可持续发展和社会安定，成为现代城市的重要安全隐患。因此，为了管理地下水资源的可持续利用和控制地面沉降，反演不同层位含水层组储水系数，了解其储水特性至关重要。

2、现有的研究估算储水系数的思路大致可以分为两类。第一类是不考虑复杂含水层系统的分层特性，将含水层组作为一个整体来研究。例如riley等用图形方法，通过绘制应力应变曲线图，根据曲线的反斜率来估计储水系数；胡谢等将insar数据和地下水位数据相结合，用谐波序列方法，通过线性回归反演储水系数；白林等使用多通道奇异谱分析方法，分离季节性变形和季节性水头变化，考虑滞后用最小二乘法反演储水系数。然而不同层位的开采情况不一样，含水层组的压缩情况不同，储水系数也不同。第二类是考虑复杂含水层系统的分层特性，利用分层标数据结合水位数据估算不同层位含水层组的储水系数。例如李江涛等依据分层标，用快速独立分量分析(fast-ica)与可变预固结水头分解方法相结合的方法，估算了分层含水层的储水系数。雷坤超等利用分层标监测站，根据预固结水头判断并估算北京平原区不同压缩层组储水系数。这类研究需要分层标数据，其实验成本高，监测点位少，仅仅能反映个别点位的储水系数情况，难以大范围推广。

3、传统的大地测量技术(如gps、水准测量等)空间分辨率低、覆盖范围小。多时相合成孔径雷达干涉测量(multi-temporal synthetic aperture radar interferometry,mt-insar)技术可以同时获取地表对微波的反射强度和相位信息，与传统测量技术相比，具有全天时、全天候、测量范围广等特点，可以在大规模尺度上监测地表变形。mt-insar技术利用同一地区的多景sar影像对时序稳定点(persistent scatterer,ps)进行精确分析，极大地降低了大气延迟等带来的测量误差，使得形变监测精度达到了厘米级到毫米级，在角反射器的辅助下甚至可实现亚毫米级监测精度，相对于传统测量方法，可提取城市区域大范围、高精度的地表三维信息和形变信息。mt-insar发展到现在比较有代表性的是以单一影像为主影像的永久散射体干涉(persistent scatterer insar，ps-insar)方法和以多幅影像为主影像的小基线(small baseline subsets,sbas)方法。

4、在空间统计分析方法中，mgwr模型描述了自变量与因变量之间的关系，还考虑了数据的空间异质性。将数据的地理位置考虑到回归参数之中，能够体现变量随空间位置的不同而出现差异。允许每个自变量拥有各自不同的空间平滑水平，解决了gwr模型所有变量同一平滑水平的缺陷，降低了估计的偏误。通过mgwr模型可以量化分层含水层地下水位变化与地面沉降的相关性，进而获得各层含水层对地面沉降的贡献率。

5、综上所述，以往估算储水系数的方法，大多在少数点位进行，或者将含水层组看作整体计算，将结果视作整个含水层组共用的储水系数。随着全球地面沉降的发育，研究出一种能够在地面沉降区，反演不同层位含水层组区域尺度储水系数，以了解其储水特性的方法至关重要。

技术实现思路

1、为解决上述问题，本发明的目的在于提供一种联合mt-insar技术和mgwr模型的分层含水层组储水系数反演方法，适用于空间尺度且不同层位的储水系数反演。该方法联合mt-insar技术和mgwr模型，量化分层含水层水位变化与地面沉降相关性，获得各层含水层对地面沉降的贡献率，根据分层含水层水位变化对地面沉降的贡献度，在空间尺度反演得到不同层位含水层组的储水系数。

2、为达到上述目的，本发明提供了一种联合mt-insar技术和mgwr模型的分层含水层组储水系数反演方法，具体包括以下步骤：

3、步骤s1：利用mt-insar处理sar影像，获取地面沉降信息，所述地面沉降信息包括地面沉降数据和地下水位数据；

4、步骤s2：利用arcgis软件的克里金插值法将需要输入mgwr模型的沉降数据和地下水位数据进行预处理；

5、步骤s3：构建mgwr模型，并将预处理后的数据输入模型，获取分层含水层地下水位变化与地面沉降的相关性；

6、步骤s4：基于分层含水层地下水位变化与地面沉降的相关性，计算分层含水层对地面沉降的贡献率；

7、步骤s5：根据得到的贡献率将mt-insar获取的地面沉降信息分层，具体为将步骤s1中获取的地面沉降信息分别乘步骤s4计算得到的分层含水层对地面沉降的贡献率，得到不同层位含水层组的地面沉降信息；

8、步骤s6：将分层的地面沉降信息结合各层位地下水位变化数据在空间尺度反演得到不同层位含水层组的储水系数。

9、在本发明一实施例中，其中，步骤s1具体包括：

10、步骤s101：设sar影像共有n+1幅，选取其中一幅作为主影像，将剩余的n幅影像作为辅影像与主影像进行配准；

11、步骤s102：利用外部获取的dem数据模拟地形相位，然后将辅影像与主影像进行差分干涉以去除地形相位和平地相位；确定基线对，并采用确定的基线组合再次进行差分干涉处理；

12、步骤s103：通过相关系数阈值法或相位离差阈值法选取永久散射体；

13、步骤s104：进行相位解缠和噪声去除处理，其中，所述噪声包括大气误差及dem引入的地形误差；

14、步骤s105：去除轨道误差并生成时间序列形变信息，其中轨道误差是指垂直轨道分量；

15、其中，在步骤s102进行干涉处理后能够得到n幅干涉相位图，然后通过式(1)对干涉相位进行分解：

16、

17、其中，为干涉相位；为平地相位，其通过在读取sar数据时利用卫星精确轨道数据去除；为地形相位，其是利用从美国srtm获取的dem数据去除；为变形相位；为大气相位，其通过sarproz软件中的aps处理消除；为aps处理过程中的热噪声和配准误差，其通过线性模型消除；

18、在消除和之后，得到雷达视距方向的变形相位然后将雷达视距方向的变形相位通过式(2)转化为垂直位移dv，即得到步骤s105的垂直轨道分量：

19、dv＝dlos/cosθ    (2)

20、其中，dlos为雷达视距方向的变形相位，θ为入射角。

21、在本发明一实施例中，其中，步骤s2是利用arcgis软件的克里金插值法将需要获取到的基于点的水位数据插值成面数据，需要计算每一层地下水位数据的变化量，具体为：

22、利用arcgis软件中的创建渔网工具建立1km×1km的格网，然后将各层地下水位变化数据和沉降数据提取均值到每个格网中，以作为模型输入数据。

23、在本发明一实施例中，其中，步骤s3构建mgwr模型的具体过程包括：

24、步骤s301：采用莫兰指数i进行分析考察每个自变量与因变量是否存在空间自相关，其中自变量为各含水层地下水位变化，因变量为累计地面沉降，具体计算过程如式(3)～式(6)：

25、

26、

27、

28、v(i)＝v(i2)-v(i)2        (6)

29、式中，n为研究对象的个数，xi和xj分别表示第i个空间单元和第j个空间单元的属性值，为所有空间单元属性值的均值，wij为空间单元i和j之间的空间权重，zscore为归一化统计的阈值，e(i)为期望值自相关，v(i)为方差；

30、其中，莫兰指数i的取值范围为[-1,1]，当莫兰指数i小于0表示地理样本之间的相关性为负空间相关性，当莫兰指数i大于0表示地理样本之间的相关性为正空间相关性，当莫兰指数i接近于0时，则表明空间分布是随机的，不存在空间自相关；

31、步骤s302：在自变量和因变量都存在空间自相关的情况下，构建mgwr模型，mgwr模型的输入包括自变量和因变量，其中，mgwr模型结构方程为：

32、

33、式中，yi为位置i处的因变量值，bwj表示第j个变量系数使用的带宽，(ui，vi)为第i个采样点的坐标，xij为第j个预测变量，βbwj(ui，vi)为位置i处第j个变量的回归系数，εi为mgwr模型在位置i处的误差项，k为第k个变量。

34、在本发明一实施例中，其中，步骤s4具体为：

35、步骤s401：根据通过mgwr模型获取的分层含水层地下水位变化与地面沉降的相关性，在沉降区域内选取分层含水层地下水位变化与地面沉降的相关系数为正的区域；

36、步骤s402：通过式(8)估算分层含水层对地面沉降的贡献率：

37、

38、式中，i为第i层含水层，ri为第i层含水层组对地面沉降的贡献率，n为含水层的数量，xi为沉降区域网格内所有正回归系数之和。

39、在本发明一实施例中，其中，步骤s6具体包括：

40、假设省略了水的可压缩性，根据含水层系统的压缩与水头变化之间的关系来计算与含水层系统可压缩性相关的骨架的蓄水量sk为：

41、

42、式中，δb为步骤s5得到不同层位含水层组的地面沉降信息，即含水层系统的变形，δh为不同层位地下水位变化数据，即水头的变化，并将估算得到的sk作为分层含水层的储水系数。

43、本发明提供联合mt-insar技术和mgwr模型的分层含水层组储水系数反演方法，与现有技术相比，至少具备以下优点：

44、1)在沉降数据获取方面，mt-insar技术可以获得大范围沉降信息，解决了现有技术采用分层标数据实验成本高，监测点位少的问题；

45、2)本发明联合mt-insar技术和mgwr模型，结合与地下水位数据的相关性，能够低成本获得区域尺度分层的沉降信息，因此可以获得分层含水层组的储水系数，极大程度的减少了所需成本，具有很强的可推广性。