CN116861647A - 联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法 - Google Patents
联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法 Download PDFInfo
- Publication number
- CN116861647A CN116861647A CN202310774616.3A CN202310774616A CN116861647A CN 116861647 A CN116861647 A CN 116861647A CN 202310774616 A CN202310774616 A CN 202310774616A CN 116861647 A CN116861647 A CN 116861647A
- Authority
- CN
- China
- Prior art keywords
- water
- mgwr
- ground
- model
- layered
- 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 134
- 238000000034 method Methods 0.000 title claims abstract description 74
- 238000005516 engineering process Methods 0.000 title claims abstract description 34
- 230000008859 change Effects 0.000 claims abstract description 42
- 239000003673 groundwater Substances 0.000 claims abstract description 40
- 238000004062 sedimentation Methods 0.000 claims abstract description 18
- 238000012545 processing Methods 0.000 claims abstract description 12
- 238000007781 pre-processing Methods 0.000 claims abstract description 5
- 230000008569 process Effects 0.000 claims description 14
- 230000001419 dependent effect Effects 0.000 claims description 13
- 238000005070 sampling Methods 0.000 claims description 7
- 230000006835 compression Effects 0.000 claims description 5
- 238000007906 compression Methods 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000006073 displacement reaction Methods 0.000 claims description 4
- 239000006185 dispersion Substances 0.000 claims description 3
- 230000008030 elimination Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 claims description 3
- 238000012544 monitoring process Methods 0.000 description 8
- 238000005305 interferometry Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 5
- 238000009499 grossing Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 230000002085 persistent effect Effects 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000012876 topography Methods 0.000 description 3
- 238000007596 consolidation process Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000000691 measurement method Methods 0.000 description 2
- 230000001932 seasonal effect Effects 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 238000009933 burial Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000001739 density measurement Methods 0.000 description 1
- 238000012880 independent component analysis Methods 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 230000006641 stabilisation Effects 0.000 description 1
- 238000011105 stabilization Methods 0.000 description 1
- 238000006467 substitution reaction 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
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)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开一种联合MT‑InSAR技术和MGWR模型的分层含水层组储水系数反演方法，其包括：S1：利用MT‑InSAR处理SAR影像，获取地面沉降信息；S2：利用ArcGIS软件的克里金插值法将需要输入MGWR模型的沉降数据和地下水位数据进行预处理；S3：构建MGWR模型，并将预处理后的数据输入模型，获取分层含水层地下水位变化与地面沉降的相关性；S4：基于分层含水层地下水位变化与地面沉降的相关性，计算分层含水层对地面沉降的贡献率；S5：根据得到的贡献率将MT‑InSAR获取的地面沉降信息分层；S6：将分层的地面沉降信息结合各层位地下水位变化数据在空间尺度反演得到不同层位含水层组的储水系数。
Description
技术领域
本发明涉及一种含水层组储水系数反演的方法，具体而言，涉及一种联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，特别是针对主要因为地下水开采导致的地面沉降区，联合MT-InSAR技术和空间统计分析方法，量化分层含水层水位变化与地面沉降相关性，根据分层含水层水位变化对地面沉降的贡献度，将沉降数据进行分层，并结合地下水位监测数据，在空间尺度反演不同层位含水层组的储水系数。
背景技术
地面沉降是指在一定的地表面积内所发生的地表海拔标高降低的现象，是一种缓变性的地质灾害，主要是由于过量开采地下水，导致地下水位降低，从而产生地面沉降的现象，造成含水层储水能力损失。不均匀的地面沉降会影响和制约当地经济建设可持续发展和社会安定，成为现代城市的重要安全隐患。因此，为了管理地下水资源的可持续利用和控制地面沉降，反演不同层位含水层组储水系数，了解其储水特性至关重要。
现有的研究估算储水系数的思路大致可以分为两类。第一类是不考虑复杂含水层系统的分层特性，将含水层组作为一个整体来研究。例如Riley等用图形方法，通过绘制应力应变曲线图，根据曲线的反斜率来估计储水系数；胡谢等将InSAR数据和地下水位数据相结合，用谐波序列方法，通过线性回归反演储水系数；白林等使用多通道奇异谱分析方法，分离季节性变形和季节性水头变化，考虑滞后用最小二乘法反演储水系数。然而不同层位的开采情况不一样，含水层组的压缩情况不同，储水系数也不同。第二类是考虑复杂含水层系统的分层特性，利用分层标数据结合水位数据估算不同层位含水层组的储水系数。例如李江涛等依据分层标，用快速独立分量分析(Fast-ICA)与可变预固结水头分解方法相结合的方法，估算了分层含水层的储水系数。雷坤超等利用分层标监测站，根据预固结水头判断并估算北京平原区不同压缩层组储水系数。这类研究需要分层标数据，其实验成本高，监测点位少，仅仅能反映个别点位的储水系数情况，难以大范围推广。
传统的大地测量技术(如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)方法。
在空间统计分析方法中，MGWR模型描述了自变量与因变量之间的关系，还考虑了数据的空间异质性。将数据的地理位置考虑到回归参数之中，能够体现变量随空间位置的不同而出现差异。允许每个自变量拥有各自不同的空间平滑水平，解决了GWR模型所有变量同一平滑水平的缺陷，降低了估计的偏误。通过MGWR模型可以量化分层含水层地下水位变化与地面沉降的相关性，进而获得各层含水层对地面沉降的贡献率。
综上所述，以往估算储水系数的方法，大多在少数点位进行，或者将含水层组看作整体计算，将结果视作整个含水层组共用的储水系数。随着全球地面沉降的发育，研究出一种能够在地面沉降区，反演不同层位含水层组区域尺度储水系数，以了解其储水特性的方法至关重要。
发明内容
为解决上述问题，本发明的目的在于提供一种联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，适用于空间尺度且不同层位的储水系数反演。该方法联合MT-InSAR技术和MGWR模型，量化分层含水层水位变化与地面沉降相关性，获得各层含水层对地面沉降的贡献率，根据分层含水层水位变化对地面沉降的贡献度，在空间尺度反演得到不同层位含水层组的储水系数。
为达到上述目的，本发明提供了一种联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，具体包括以下步骤：
步骤S1：利用MT-InSAR处理SAR影像，获取地面沉降信息，所述地面沉降信息包括地面沉降数据和地下水位数据；
步骤S2：利用ArcGIS软件的克里金插值法将需要输入MGWR模型的沉降数据和地下水位数据进行预处理；
步骤S3：构建MGWR模型，并将预处理后的数据输入模型，获取分层含水层地下水位变化与地面沉降的相关性；
步骤S4：基于分层含水层地下水位变化与地面沉降的相关性，计算分层含水层对地面沉降的贡献率；
步骤S5：根据得到的贡献率将MT-InSAR获取的地面沉降信息分层，具体为将步骤S1中获取的地面沉降信息分别乘步骤S4计算得到的分层含水层对地面沉降的贡献率，得到不同层位含水层组的地面沉降信息；
步骤S6：将分层的地面沉降信息结合各层位地下水位变化数据在空间尺度反演得到不同层位含水层组的储水系数。
在本发明一实施例中，其中，步骤S1具体包括：
步骤S101：设SAR影像共有N+1幅，选取其中一幅作为主影像，将剩余的N幅影像作为辅影像与主影像进行配准；
步骤S102：利用外部获取的DEM数据模拟地形相位，然后将辅影像与主影像进行差分干涉以去除地形相位和平地相位；确定基线对，并采用确定的基线组合再次进行差分干涉处理；
步骤S103：通过相关系数阈值法或相位离差阈值法选取永久散射体；
步骤S104：进行相位解缠和噪声去除处理，其中，所述噪声包括大气误差及DEM引入的地形误差；
步骤S105：去除轨道误差并生成时间序列形变信息，其中轨道误差是指垂直轨道分量；
其中，在步骤S102进行干涉处理后能够得到N幅干涉相位图，然后通过式(1)对干涉相位进行分解：
其中，为干涉相位；/>为平地相位，其通过在读取SAR数据时利用卫星精确轨道数据去除；/>为地形相位，其是利用从美国SRTM获取的DEM数据去除；/>为变形相位；/>为大气相位，其通过SARPROZ软件中的APS处理消除；/>为APS处理过程中的热噪声和配准误差，其通过线性模型消除；
在消除和/>之后，得到雷达视距方向的变形相位/>然后将雷达视距方向的变形相位通过式(2)转化为垂直位移dv，即得到步骤S105的垂直轨道分量：
dv＝dlos/cosθ (2)
其中，dlos为雷达视距方向的变形相位，θ为入射角。
在本发明一实施例中，其中，步骤S2是利用ArcGIS软件的克里金插值法将需要获取到的基于点的水位数据插值成面数据，需要计算每一层地下水位数据的变化量，具体为：
利用ArcGIS软件中的创建渔网工具建立1km×1km的格网，然后将各层地下水位变化数据和沉降数据提取均值到每个格网中，以作为模型输入数据。
在本发明一实施例中，其中，步骤S3构建MGWR模型的具体过程包括：
步骤S301：采用莫兰指数I进行分析考察每个自变量与因变量是否存在空间自相关，其中自变量为各含水层地下水位变化，因变量为累计地面沉降，具体计算过程如式(3)～式(6)：
V(I)＝V(I2)-V(I)2 (6)
式中，n为研究对象的个数，xi和xj分别表示第i个空间单元和第j个空间单元的属性值，为所有空间单元属性值的均值，wij为空间单元i和j之间的空间权重，zscore为归一化统计的阈值，E(I)为期望值自相关，V(I)为方差；
其中，莫兰指数I的取值范围为[-1,1]，当莫兰指数I小于0表示地理样本之间的相关性为负空间相关性，当莫兰指数I大于0表示地理样本之间的相关性为正空间相关性，当莫兰指数I接近于0时，则表明空间分布是随机的，不存在空间自相关；
步骤S302：在自变量和因变量都存在空间自相关的情况下，构建MGWR模型，MGWR模型的输入包括自变量和因变量，其中，MGWR模型结构方程为：
式中，yi为位置i处的因变量值，bwj表示第j个变量系数使用的带宽，(ui，vi)为第i个采样点的坐标，xij为第j个预测变量，βbwj(ui，vi)为位置i处第j个变量的回归系数，εi为MGWR模型在位置i处的误差项，k为第k个变量。
在本发明一实施例中，其中，步骤S4具体为：
步骤S401：根据通过MGWR模型获取的分层含水层地下水位变化与地面沉降的相关性，在沉降区域内选取分层含水层地下水位变化与地面沉降的相关系数为正的区域；
步骤S402：通过式(8)估算分层含水层对地面沉降的贡献率：
式中，i为第i层含水层，Ri为第i层含水层组对地面沉降的贡献率，n为含水层的数量，xi为沉降区域网格内所有正回归系数之和。
在本发明一实施例中，其中，步骤S6具体包括：
假设省略了水的可压缩性，根据含水层系统的压缩与水头变化之间的关系来计算与含水层系统可压缩性相关的骨架的蓄水量SK为：
式中，Δb为步骤S5得到不同层位含水层组的地面沉降信息，即含水层系统的变形，Δh为不同层位地下水位变化数据，即水头的变化，并将估算得到的SK作为分层含水层的储水系数。
本发明提供联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，与现有技术相比，至少具备以下优点：
1)在沉降数据获取方面，MT-InSAR技术可以获得大范围沉降信息，解决了现有技术采用分层标数据实验成本高，监测点位少的问题；
2)本发明联合MT-InSAR技术和MGWR模型，结合与地下水位数据的相关性，能够低成本获得区域尺度分层的沉降信息，因此可以获得分层含水层组的储水系数，极大程度的减少了所需成本，具有很强的可推广性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明一实施例中流程示意图；
图2为本发明一实施例中的MT-InSAR技术流程示意图；
图3为本发明一实施例中的MGWR模型不同自变量空间平滑水平原理示意图；
图4为本发明一实施例中MGWR模型运行步骤示意图；
图5为本发明一实施例获取分层含水层对地面沉降贡献率流程示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为了语言的简洁及表述清晰，在本发明说明书中，相关术语可能会以缩略语的形式出现，现将本发明涉及的缩略语的中英文全称说明如下，以利完全了解：
MT-InSAR：永久散射体合成孔径雷达干涉测量，Multi-temporal Syntheticaperture radar interferometry；
SAR：合成孔径雷达，Synthetic Aperture Radar；
SLC：单视复数影像，Single Look Complex；
SRTM：航天飞机雷达地形测绘使命，Shuttle Radar Topography Mission；
DEM：数字高程模型，Digital Elevation Model；
LOS：视线，Line of Sight；
APS：大气相位分析，Atmospheric phase analysis；
MGWR：多尺度地理加权回归，Multiscale Geographically WeightedRegression；
AICc：赤池信息量准则，Akaike Information Criterions；
本发明所欲解决的现有技术的不足主要包括：
1)较低成本获得区域尺度的沉降量。以往沉降数据是通过分层标数据获取，需要的分层标数据，实验成本高，监测点位少，仅仅能反映个别点位的储水系数情况，难以大范围推广。本发明通过MT-InSAR技术可以获得大范围沉降信息，获得区域尺度的沉降量。
2)通过贡献率获得分层的储水系数。以往不考虑复杂含水层系统的分层特性，将含水层组作为一个整体来研究。本发明通过MGWR模型可以量化不同层含水层组与地面沉降相关性，根据不同层位含水层组对地面沉降的贡献率，获得分层的沉降量，因此可以获得分层含水层组的储水系数。
本发明提供的联合MT-InSAR技术和MGWR模型的分层含水层组区域尺度储水系数反演方法，能够获取空间尺度不同层位含水层组储水系数。其中MGWR模型需要的输入数据有：地面沉降值、同期不同层位地下水位变化数据。本发明的主要技术路线为：先对MT-InSAR技术获取的地面沉降数据和地下水位监测数据进行预处理，获取MGWR模型所需的输入数据；数据输入模型后获取不同层位含水层水位变化与地面沉降的相关性，进一步求解分层含水层水位变化对地面沉降的贡献度；最后根据贡献度将MT-InSAR技术获取的地面沉降数据分层到各含水层组，进而获取分层含水层组区域尺度储水系数。以下将通过具体实施方式进行详细说明。
本实施例提供一种联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，具体包括以下步骤：
步骤S1：利用MT-InSAR处理SAR影像，获取地面沉降信息，所述地面沉降信息包括地面沉降数据和地下水位数据(同期不同层位地下水位变化数据)；
MT-InSAR是一种利用多幅SAR影像进行高精度、高密度测量点提取分析的技术。InSAR(合成孔径雷达干涉测量)可以获取大范围的地表高程信息，而通过差分InSAR技术(D-InSAR)可以得到地表微小形变。在实际应用中，形变的发生是缓慢而连续的，因此D-InSAR技术只能获取某一时刻到另外一个时刻的形变，却无法得知这一形变过程的发展演化，并且若时间间隔太长则容易引起时间失相干，进而无法获得形变。为此，Ferriti等在1999年提出了永久散射体InSAR测量技术(persistent scatterer InSAR,PS-InSAR)，即最初的MT-InSAR方法。这种方法首先确定地表永久散射体，并对其进行持续的干涉测量，以获取对应位置的高程信息，随后通过相位分析，分离出目标点的形变时序结果。小基线集InSAR技术(SBAS-InSAR)是另外一种MT-InSAR方法，其与PS-InSAR最大的区别就是采用多主影像的方式完成差分干涉处理。
图2为本发明一实施例中的MT-InSAR技术流程示意图，如图2所示，在本实施例中，其中，步骤S1具体包括：
步骤S101：设SAR影像共有N+1幅，选取其中一幅(通常选择中间时刻获取的影像)作为主影像，将剩余的N幅影像作为辅影像与主影像进行配准；
步骤S102：利用外部获取的DEM数据模拟地形相位，然后将辅影像与主影像进行差分干涉以去除地形相位和平地相位；确定基线对，并采用确定的基线组合再次进行差分干涉处理；
在本实施例中，若采用基于PS-InSAR的MT-InSAR方法，则基线对为单个主影像与多个辅影像组成配对；若采用基于SBAS-InSAR的MT-InSAR方法，则采用类似排列组合的方式，在保证相干性的情况下允许每个影像都参与配对。
步骤S103：通过相关系数阈值法或相位离差阈值法选取永久散射体(persistentscatterer，PS)，即建立PS点形变反演模型；
步骤S104：进行相位解缠和噪声去除处理，其中，噪声通常包括大气误差、DEM引入的地形误差等；
步骤S105：去除轨道误差并生成时间序列形变信息，其中轨道误差是指垂直轨道分量；在实际执行时后期的误差去除，需要根据实际情况来识图判别，有可能要进行多次的误差去除；
其中，在步骤S102进行干涉处理后能够得到N幅干涉相位图，然后通过式(1)对干涉相位进行分解：
其中，为干涉相位；/>为平地相位，其通过在读取SAR数据时利用卫星精确轨道数据去除；/>为由地形起伏所造成的地形相位，其是利用从美国SRTM(ShuttleRadar Topography Mission)获取的DEM数据去除；/>为两次采集图像时由于地面位移而引起的变形相位；/>为由于大气组分的贡献而产生的大气相位，其可通过SARPROZ软件中的APS(Atmospheric phase analysis)处理消除；/>为APS处理过程中的热噪声和配准误差，其可通过线性模型消除；
在消除和/>之后，可以得到雷达视距LOS(Line ofSight)方向(包括水平方向和垂直方向)的变形相位/>然后将雷达视距方向的变形相位通过式(2)转化为垂直位移dv，即得到步骤S105的垂直轨道分量：
dv＝dlos/cosθ (2)
其中，dlos为雷达视距方向的变形相位，θ为入射角。
步骤S2：利用ArcGIS(一种GIS平台)软件的克里金插值法将需要输入MGWR模型的沉降数据和地下水位数据进行预处理；
在本实施例中，其中，由于获取到的地下水位数据通常是点要素或者是等值线要素，因此步骤S2利用ArcGIS软件的克里金插值法将需要获取到的基于点的水位数据插值成面数据，由于地质结构的差异性，地下水位的埋深与地面沉降并不直接具有相关性，因此需要计算每一层地下水位数据的变化量，具体为：
利用ArcGIS软件中的创建渔网工具建立1km×1km的格网，然后将各层地下水位变化数据和沉降数据提取均值到每个格网中，以作为模型输入数据，这种方式可以减少数据的冗余。
步骤S3：构建MGWR模型，并将预处理后的数据输入模型，获取分层含水层地下水位变化与地面沉降的相关性；
在本实施例中，其中，步骤S3构建MGWR模型的具体过程包括：
步骤S301：采用莫兰指数(Moran’s)I进行分析考察每个自变量与因变量是否存在空间自相关，其中自变量为各含水层地下水位变化，因变量为累计地面沉降，具体计算过程如式(3)～式(6)：
V(I)＝V(I2)-V(I)2 (6)
式中，n为研究对象的个数，xi和xj分别表示第i个空间单元和第j个空间单元的属性值，为所有空间单元属性值的均值，wij为空间单元i和j之间的空间权重，zscore为归一化统计的阈值，E(I)为期望值自相关，V(I)为方差；
其中，莫兰指数I的取值范围为[-1,1]，当莫兰指数I小于0表示地理样本之间的相关性为负空间相关性，当莫兰指数I大于0表示地理样本之间的相关性为正空间相关性，当莫兰指数I接近于0时，则表明空间分布是随机的，不存在空间自相关；
步骤S302：在自变量和因变量都存在空间自相关的情况下，构建MGWR(MultiscaleGeographic Weighted Regression)模型，MGWR模型的输入包括自变量和因变量，图3为本发明一实施例中的MGWR模型不同自变量空间平滑水平原理示意图，如图3所示，MGWR允许每个自变量拥有各自不同的空间平滑水平，以降低估计的偏误，其中，MGWR模型结构方程为：
式中，yi为位置i处的因变量值，bwj表示第j个变量系数使用的带宽，(ui，vi)为第i个采样点的坐标，xij为第j个预测变量，βbwj(ui，vi)为位置i处第j个变量的回归系数，εi为MGWR模型在位置i处的误差项，k为第k个变量。
其中，核函数和带宽选择例如可以分别使用高斯权函数和AICc(AkaikeInformation Criterions，赤池信息量)准则，但并不限于此，在其他实施例可以根据需求选择。
图4为本发明一实施例中MGWR模型运行步骤示意图，如图4所示，以四层含水层地下水位变化量为自变量，累计地面沉降为因变量，通过MGWR模型进行计算，先确定一因变量的带宽，然后通过AICc准则确定样点周围样点的权重举证，之后根据带宽确定核函数，以得到每个样点的回归方程，再重复确定每个因变量的带宽以重复计算每个样点的回归方程，进而得到地面沉降与各层含水层组水位变化的回归系数。
步骤S4：基于分层含水层地下水位变化与地面沉降的相关性，计算分层含水层对地面沉降的贡献率；
图5为本发明一实施例获取分层含水层对地面沉降贡献率流程示意图，如图5所示，在本实施例中，其中，步骤S4具体为：
步骤S401：根据通过MGWR模型获取的分层含水层地下水位变化与地面沉降的相关性，在沉降区域内选取分层含水层地下水位变化与地面沉降的相关系数为正的区域，即正向影响的区域；
步骤S402：通过式(8)估算分层含水层对地面沉降的贡献率：
式中，i为第i层含水层，Ri为第i层含水层组对地面沉降的贡献率，n为含水层的数量，xi为沉降区域网格内所有正回归系数之和。
步骤S5：根据得到的贡献率将MT-InSAR获取分层的地面沉降信息，具体为将步骤S1中获取的地面沉降信息分别乘步骤S4计算得到的分层含水层对地面沉降的贡献率，得到不同层位含水层组的地面沉降信息；
步骤S6：将分层的地面沉降信息结合各层位地下水位变化数据在空间尺度反演得到不同层位含水层组的储水系数。
承压含水层的特定储存系数S是指单位水头下降时每单位含水层体积从可压缩含水层中排出的水量，通常，含水层系统的可压缩性远大于水。因此在本实施例中，其中，步骤S6具体包括：
假设省略了水的可压缩性，根据含水层系统的压缩与水头变化之间的关系来计算与含水层系统可压缩性相关的骨架的蓄水量SK为：
式中，Δb为步骤S5得到不同层位含水层组的地面沉降信息，即含水层系统的变形，Δh为不同层位地下水位变化数据，即水头的变化，并将估算得到的SK作为分层含水层的储水系数。
本发明提供联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，与现有技术相比，至少具备以下优点：
1)在沉降数据获取方面，MT-InSAR技术可以获得大范围沉降信息，解决了现有技术采用分层标数据实验成本高，监测点位少的问题；
2)本发明联合MT-InSAR技术和MGWR模型，结合与地下水位数据的相关性，能够低成本获得区域尺度分层的沉降信息，因此可以获得分层含水层组的储水系数，极大程度的减少了所需成本，具有很强的可推广性。
本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。
本领域普通技术人员可以理解：实施例中的装置中的模块可以按照实施例描述分布于实施例的装置中，也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。
Claims (6)
1.一种联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，其特征在于，包括以下步骤：
步骤S1：利用MT-InSAR处理SAR影像，获取地面沉降信息，所述地面沉降信息包括地面沉降数据和地下水位数据；
步骤S2：利用ArcGIS软件的克里金插值法将需要输入MGWR模型的沉降数据和地下水位数据进行预处理；
步骤S3：构建MGWR模型，并将预处理后的数据输入模型，获取分层含水层地下水位变化与地面沉降的相关性；
步骤S4：基于分层含水层地下水位变化与地面沉降的相关性，计算分层含水层对地面沉降的贡献率；
步骤S5：根据得到的贡献率将MT-InSAR获取的地面沉降信息分层，具体为将步骤S1中获取的地面沉降信息分别乘步骤S4计算得到的分层含水层对地面沉降的贡献率，得到不同层位含水层组的地面沉降信息；
步骤S6：将分层的地面沉降信息结合各层位地下水位变化数据在空间尺度反演得到不同层位含水层组的储水系数。
2.根据权利要求1所述的联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，其特征在于，步骤S1具体包括：
步骤S101：设SAR影像共有N+1幅，选取其中一幅作为主影像，将剩余的N幅影像作为辅影像与主影像进行配准；
步骤S102：利用外部获取的DEM数据模拟地形相位，然后将辅影像与主影像进行差分干涉以去除地形相位和平地相位；确定基线对，并采用确定的基线组合再次进行差分干涉处理；
步骤S103：通过相关系数阈值法或相位离差阈值法选取永久散射体；
步骤S104：进行相位解缠和噪声去除处理，其中，所述噪声包括大气误差及DEM引入的地形误差；
步骤S105：去除轨道误差并生成时间序列形变信息，其中轨道误差是指垂直轨道分量；
其中，在步骤S102进行干涉处理后能够得到N幅干涉相位图，然后通过式(1)对干涉相位进行分解：
其中，为干涉相位；/>为平地相位，其通过在读取SAR数据时利用卫星精确轨道数据去除；/>为地形相位，其是利用从美国SRTM获取的DEM数据去除；/>为变形相位；为大气相位，其通过SARPROZ软件中的APS处理消除；/>为APS处理过程中的热噪声和配准误差，其通过线性模型消除；
在消除和/>之后，得到雷达视距方向的变形相位/>然后将雷达视距方向的变形相位通过式(2)转化为垂直位移dv，即得到步骤S105的垂直轨道分量：
dv＝dlos/cosθ (2)
其中，dlos为雷达视距方向的变形相位，θ为入射角。
3.根据权利要求1所述的联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，其特征在于，步骤S2是利用ArcGIS软件的克里金插值法将需要获取到的基于点的水位数据插值成面数据，需要计算每一层地下水位数据的变化量，具体为：
利用ArcGIS软件中的创建渔网工具建立1km×1km的格网，然后将各层地下水位变化数据和沉降数据提取均值到每个格网中，以作为模型输入数据。
4.根据权利要求1所述的联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，其特征在于，步骤S3构建MGWR模型的具体过程包括：
步骤S301：采用莫兰指数I进行分析考察每个自变量与因变量是否存在空间自相关，其中自变量为各含水层地下水位变化，因变量为累计地面沉降，具体计算过程如式(3)～式(6)：
V(I)＝V(I2)-V(I)2 (6)
式中，n为研究对象的个数，xi和xj分别表示第i个空间单元和第j个空间单元的属性值，为所有空间单元属性值的均值，wij为空间单元i和j之间的空间权重，zscore为归一化统计的阈值，E(I)为期望值自相关，V(I)为方差；
其中，莫兰指数I的取值范围为[-1,1]，当莫兰指数I小于0表示地理样本之间的相关性为负空间相关性，当莫兰指数I大于0表示地理样本之间的相关性为正空间相关性，当莫兰指数I接近于0时，则表明空间分布是随机的，不存在空间自相关；
步骤S302：在自变量和因变量都存在空间自相关的情况下，构建MGWR模型，MGWR模型的输入包括自变量和因变量，其中，MGWR模型结构方程为：
式中，yi为位置i处的因变量值，bwj表示第j个变量系数使用的带宽，(ui，vi)为第i个采样点的坐标，xij为第j个预测变量，βbwj(ui，vi)为位置i处第j个变量的回归系数，εi为MGWR模型在位置i处的误差项，k为第k个变量。
5.根据权利要求1所述的联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，其特征在于，步骤S4具体为：
步骤S401：根据通过MGWR模型获取的分层含水层地下水位变化与地面沉降的相关性，在沉降区域内选取分层含水层地下水位变化与地面沉降的相关系数为正的区域；
步骤S402：通过式(8)估算分层含水层对地面沉降的贡献率：
式中，i为第i层含水层，Ri为第i层含水层组对地面沉降的贡献率，n为含水层的数量，xi为沉降区域网格内所有正回归系数之和。
6.根据权利要求1所述的联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法，其特征在于，步骤S6具体包括：
假设省略了水的可压缩性，根据含水层系统的压缩与水头变化之间的关系来计算与含水层系统可压缩性相关的骨架的蓄水量SK为：
式中，Δb为步骤S5得到不同层位含水层组的地面沉降信息，即含水层系统的变形，Δh为不同层位地下水位变化数据，即水头的变化，并将估算得到的SK作为分层含水层的储水系数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310774616.3A CN116861647A (zh) | 2023-06-28 | 2023-06-28 | 联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310774616.3A CN116861647A (zh) | 2023-06-28 | 2023-06-28 | 联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116861647A true CN116861647A (zh) | 2023-10-10 |
Family
ID=88224405
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310774616.3A Pending CN116861647A (zh) | 2023-06-28 | 2023-06-28 | 联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116861647A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117593542A (zh) * | 2023-11-27 | 2024-02-23 | 首都师范大学 | 粮食生产区地面沉降差异演化特征计算方法、装置及介质 |
-
2023
- 2023-06-28 CN CN202310774616.3A patent/CN116861647A/zh active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117593542A (zh) * | 2023-11-27 | 2024-02-23 | 首都师范大学 | 粮食生产区地面沉降差异演化特征计算方法、装置及介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111142119B (zh) | 一种基于多源遥感数据的矿山地质灾害动态识别与监测方法 | |
| Brasington et al. | Monitoring and modelling morphological change in a braided gravel‐bed river using high resolution GPS‐based survey | |
| Wöppelmann et al. | Is land subsidence increasing the exposure to sea level rise in Alexandria, Egypt? | |
| CN104931022A (zh) | 基于星载激光测高数据的卫星影像立体区域网平差方法 | |
| CN109212522B (zh) | 一种基于双基星载sar的高精度dem反演方法及装置 | |
| Mason et al. | Measurement of recent intertidal sediment transport in Morecambe Bay using the waterline method | |
| Goetz et al. | Quantifying uncertainties in snow depth mapping from structure from motion photogrammetry in an alpine area | |
| CN113960595A (zh) | 一种地表形变监测方法及系统 | |
| Li et al. | Intertidal topographic maps and morphological changes in the German Wadden Sea between 1996–1999 and 2006–2009 from the waterline method and SAR images | |
| CN116861647A (zh) | 联合MT-InSAR技术和MGWR模型的分层含水层组储水系数反演方法 | |
| Yim et al. | Tsunami modeling, fluid load simulation, and validation using geospatial field data | |
| CN111257870B (zh) | 一种利用InSAR监测数据的采煤沉陷积水区水下地形反演方法 | |
| CN113238228B (zh) | 基于水准约束的三维地表形变获取方法、系统及装置 | |
| Feng et al. | A hierarchical network densification approach for reconstruction of historical ice velocity fields in East Antarctica | |
| CN113900069A (zh) | 一种基于干涉成像高度计的垂线偏差计算方法及其系统 | |
| Wang et al. | Applicability of an ultra-long-range terrestrial laser scanner to monitor the mass balance of Muz Taw Glacier, Sawir Mountains, China | |
| Ardha et al. | Utilization of Sentinel-1 satellite imagery data to support land subsidence analysis in DKI Jakarta, Indonesia. | |
| Javadnejad et al. | Unmanned aircraft systems-based photogrammetry for ground movement monitoring | |
| CN114046774B (zh) | 综合cors网和多源数据的地面形变连续监测方法 | |
| Aidi et al. | UAV-based gully retrogressive erosion status dynamic variability investigations in Chinese Loess Plateau | |
| Hu et al. | Time-series spaceborne InSAR monitoring of permafrost deformation combined with backscattering characteristics: Case study from the Qinghai–Tibet engineering corridor | |
| Gesch | Topography-based analysis of Hurricane Katrina inundation of New Orleans | |
| Li et al. | Detecting, Monitoring, and Analyzing the Surface Subsidence in the Yellow River Delta (China) Combined with CenterNet Network and SBAS‐InSAR | |
| Eskandari et al. | Validation of Full-Resolution Dinsar-Derived Vertical Displacement in Cultural Heritage Monitoring: Integration with Geodetic Levelling Measurements | |
| DE MARCO | Monitoring and modelling cryosphere processes using high resolution surveys |
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 |