地面沉降－回弹及地下水位波动的InSAR长时序监测——以德州市为例

葛大庆，殷跃平，王 艳，张 玲，郭小方，王 毅

(1.中国地质大学(北京)水资源与环境学院，北京 100083;2.中国国土资源航空物探遥感中心，北京 100083;3.中国地质环境监测院，北京 100081)

0 引言

长期过量开采地下水，特别是深层地下水，是目前我国一些平原和三角洲地区产生地面沉降的主要原因[1]。地下水位的持续下降改变了原有的补径排特征，形成了区域性地下水位降落漏斗，地面沉降随之产生并扩大。我国已形成了以华北平原、长江三角洲及汾渭断陷盆地为主的三大地面沉降严重区，华北平原已出现近20个大型复合水位降落漏斗，地面沉降呈现出区域性且连片分布的特征[2－3]。

雷达干涉测量(InSAR)技术为区域性地面沉降的快速准确和连续监测提供了重要手段。以永久性散射体干涉测量(permanent scatterer InSAR，PSIn-SAR)[4]为代表的相干目标(coherent target，CT)时间序列分析技术克服了差分InSAR技术所面临的失相干和大气相位延迟等问题，通过对相干目标的差分干涉相位序列进行时序分析，逐个估计和分离地形相位、形变相位、大气相位和噪声等信息，进而解算相干目标的形变速率和累积形变量。与之类似的还有短基线集(small baseline subset，SBAS)[5]以及干涉测量点目标分析(interferometric point target analysis，IPTA)[6]等技术，这些技术已在地面沉降、滑坡、开采沉陷、地震和冰川滑移等多尺度缓慢地表形变研究中得以发展和应用[7－9]。

德州市位于山东省北部，是华北平原地面沉降和深层地下水位降落漏斗极为典型的地区之一。其地下水位降落漏斗在1978年前已经形成，是目前华北平原水位降落深度最大的漏斗，与衡水、沧州等漏斗连成一片。该区的地面沉降在1989年以前就已经出现，至2010年累积沉降量超过1 m。与地面沉降对应的是，德州市地下水位呈现出多年连续下降和年内季节性波动的特征[10]。由于地面沉降与地下水位下降相关，因此在德州市深、浅层地下水位均呈现出季节性波动的情况下，地面沉降场的时空演变特征值得研究，这包括:沉降场的空间变化特征，即沉降中心的移动及其变化幅度;沉降场的时间变化特征，即季节性波动，包括是否会出现地面回弹，回弹幅度多大，地面沉降与地下水位波动是否同步等方面。针对上述问题，本文以相干目标InSAR时序分析方法为主，利用2004年1月—2010年10月近7 a的ENVISAT卫星数据研究分析德州地区地面沉降－回弹的动态变化特征及其与地下水位波动的对应关系，以揭示该地区地面沉降季节性变化的产生原因和诱发机制。

1 相干目标InSAR时序分析方法

InSAR时序分析主要有PSInSAR和SBAS两类方法。根据干涉像对的组合模式，前者以单一主影像构成差分干涉图序列，后者则利用短基线原则以多个主影像构成若干个干涉纹图子集。单一主影像的干涉纹图序列处理主要顾及大气相位和非线性形变的统计模型，有利于对大气相位的准确估计，但对数据量要求过高，且受相干性影响而对大变形的监测能力不足;多个主影像则降低了对数据量的要求，将时空基线小于一定阈值的干涉像对组合，生成干涉纹图序列，增加对同一信号的采样密度，有利于准确求解形变速率，在数据量较小的情况下也可以实现。在选择相干目标的方式上，前者以点目标为主，后者兼顾点目标和分布式目标[5]。本文集成了2种方法的优点，以短基线为准则构成差分干涉相位图，利用点目标识别算法提取相干目标。以Delanay三角网连接相邻的相干像元，用二维周期图估计点间形变速率和高程误差改正。以此为基础，根据大气相位特征，利用奇异值分解算法解算空间滤波后的残余相位，并对解算结果进行时域滤波，以求解非线性形变，得到不同时刻的累积变形量，再利用线性回归法求解每个相干目标的变形速率。

1.1 干涉相位时序分析的基本模型

对于给定的M景雷达数据，在短基线条件下生成N个差分干涉图(一般情况下，N＞M)，对于其中的任一差分干涉图k(k=1，2，…，N)，其相位模型为

显然，待求解量d与InSAR参数不直接相关，但受其他分量的影响。对于大气分量，其本质上是2次大气延迟影响的叠加，在空间上具有一定的相关性[4]。大气波动在时间上的随机特性决定了2次互差也为随机量。需要重点区分的是大气相位与轨道误差，虽然二者都具有空间低频的特性，但各自的影响范围有所不同。基线误差多是全局性的，而大气影响半径有限，对于大范围处理，二者应区别对待。高程误差与垂直基线相关，垂直基线是变量。变形量与干涉图的时间基线相关，时间间隔是变量，相对变形的速率与时间相关。

1.2 地表形变参数的提取算法

1.2.1 线性分量的估计

考虑到大气的空间相关性，如果对相邻2点(小于大气相关距离)求差，则可削弱大气影响。相邻相干目标i和j差分干涉相位的互差△为

由于相邻目标间的相对变形速率与时间间隔相关，干涉图的时间基线为自变量，因而可将形变量d分解为线性形变速率和非线性形变量2部分。如此，时序分析过程则转为参数估计，即

式中:CB为与垂直基线相关的系数;T为时间基线;△ε为相对高程误差;△υ为相对形变速率;μNL为非线性形变量;α为大气相位;n为噪声。

对于缠绕相位，其周期数由相对高程误差和线性形变速率共同决定，采用二维参数估计方法可实现对相位周期的估算。由此，构建目标函数

并将其从相位互差中减去，得到残余相位，即

显然，在准确估计目标函数的参数△ε和△υ时，残余相位将最小化。

由于相邻2点存在多个干涉图，即存在N维B和T，因而，对△ε和△v的估计使模型相关系数最大化，即

在缠绕相位条件下，△ε和△v为周期函数的二维频率，可利用二维周期图估计使模型相关系数最大化;在解缠相位条件下，则转换为二维线性函数，利用参数估计方法可实现对△ε和△v的求解。

在相邻2点解算时需对所有相干目标建立连接关系，整体求解形变场的速率和高程误差改正量。这一目标的实现可通过建立Delanay三角网的方法完成，也可以利用冗余网构建更为复杂的连接关系，强化对待解算方程组的约束。利用邻近法则将所有距离满足大气相关距离的相干目标连接起来，在求解完成相邻点间的互差后，通过最小二乘或加权平均的方法求解每个目标相对于参考点的变形速率，即总体形变速率场。

1.2.2 非线性分量的估计

对于具有显著非线性形变过程，仍需对残余相位进行更为复杂的处理，以提取非线性形变量。从差分相位中去除高程误差估计值和线性速率后，残余相位为

由于已解算出整周相位，因而残余相位为相位主值，其大小

这里，首先要实现对单个差分干涉图中残余相位的滤波处理。由于大气表现出空间低频特性，而非线性变形的空间变化范围较小，相对大气而言表现为高通特性。因而，对残余相位图进行空间低通滤波可以进一步弱化大气影响。在短基线集条件下求解的非线性形变量，不同于PSInSAR经典处理策略，不直接对大气进行估计，而是通过滤波减弱大气的影响。残余相位经空间高通滤波后，大气分量已经减弱。此时，进一步的处理是求解与雷达数据对应的不同时刻的非线性变形量，包含噪声等影响。根据主辅影像的关系，可将解缠后的残余相位分解为

式中:p和q分别表示生成第k(k=1，…，N)景差分干涉图的主辅影像的获取时间。由于采用短基线原则，在求解每个时刻对应的残余相位时出现秩亏方程组，解决这一问题的办法就是奇异值分解法(singular value decomposition，SVD)[3]。在此基础上，对M景雷达数据的残余相位进行时域低通滤波处理，以提取最终的非线性形变量，可表示为

在求解出非线性形变量后，每个相干目标对应的相变序列为

式中vest为线性形变速率。

1.2.3 平均形变速率的重新估计

对于显著的非线性形变过程，利用前述二维线性模型估计出来的线性速率往往偏小。在恢复整个形变序列的基础上，可根据对应的每个目标的形变序列，利用最小二乘法重新估计形变速率，即

式中d和T分别为相对起始时刻的形变量和时间间隔。

2 雷达数据及InSAR时序分析处理

研究中共获取2004年1月—2010年10月近7 a的ENVISAT卫星ASAR数据。该数据覆盖了德州地区主要沉降区，总数据量为46景(Track－447，Frame－2853)，其基线分布如图1所示。

图1 差分干涉图序列的短基线关系Fig.1 Small baseline of the differential interferograms series

需要说明的是，利用幅度离散指数[1]、相干系数[5]以及子视相关[6]均可以很好地识别出相干目标候选点。对于候选点，只有在满足模型相关函数最大化约束条件下才能称为时序分析中的相干目标。实际上，一些目标后向散射特性稳定，但非线性变化强烈，不满足该条件。因而，可通过降低相关系数的阈值，将其纳入到形变解算网络中，以增加监测点的密度。

时空基线阈值的设置取决于SAR系统参数。空间基线主要考虑雷达系统临界基线，而时间基线则应根据变形幅度大小，短时间基线相当于降低了形变场的梯度，便于相位解缠和变形参数估计。因而，本研究中设定的短基线条件为:时间基线＜450 d，空间基线＜300 m。

围绕德州地区地面沉降和地下水位漏斗的基本分布特征，分别针对区域性大范围沉降和重点沉降中心进行了2类处理:利用2007—2010年4 a间的数据提取了德州地区100 km×100 km平均地面沉降速率;分析全区地面沉降和地下水位变化关系，并确定时序分析过程中稳定的地面参考点。在此基础上，对德州城区内20 km×20 km范围的主要沉降中心的7 a数据进行了时序分析处理，提取每个相干目标的沉降率和累积沉降量序列，用以详细分析多年和年度内沉降中心的时空变化特征。

3 区域性地面沉降与地下水位变化

3.1 区域地面沉降分布

图2示出德州地区约50 km×60 km范围内2007年初—2010年底4 a间平均地面沉降速率。

图2 德州地区2007—2010年CTInSAR监测平均地面沉降率Fig.2 Average subsidence velocity of Dezhou area from 2007 to 2010 derived by CTInSAR data

从区域分布上来看，德州市区沉降较为显著，沉降区向西北方向延伸，与河北省景县和故城相连，构成连片沉降区。区内景县沉降向西北方向扩展，是目前衡水—德州沉降区的重要组成部分，沉降速率均达到40～50 mm/a。德州沉降区向东南方向，分布着3个明显的沉降中心，分别为陵县、临邑县和平原县，其中以陵县最为突出，中心位于县城西北菜园村，最大沉降速率达75 mm/a。临邑县沉降中心位于临盘采油厂，中心处沉降速率达50～60 mm/a。

3.2 与深层地下水位漏斗的比较分析

InSAR监测结果准确揭示了德州地区区域性沉降和单个漏斗的分布。连续多年的地下水位变化监测结果表明，德州市深层地下水位降落漏斗是目前鲁北地区最大的漏斗，与周边的小漏斗相连接，引起区域性地面沉降[8－9]。深层地下水开采(第Ⅳ，Ⅴ含水层组)是该地区地面沉降的主要原因。以临邑为例，该沉降中心深层地下水降落漏斗的－30 m等水压线与德城区深层漏斗相连，已成为德州漏斗的一部分;－40 m等水压线包围临邑县城和临盘镇，圈闭面积达192 km2。漏斗中心在临盘采油厂一带，水位埋深超过59.70 m，水位标高－41.82 m左右。图2中的沉降中心准确地反映了目前该漏斗的分布位置和影响范围。由于本区地下水水位动态属于较稳定的连续开采消耗型，年水位变化趋势与开采量大小密切相关。临邑深层地下水开采量达3.56×104m3/a。其中，第Ⅳ和第Ⅴ含水层组地下水开采量约占90%，地下水位多年平均降速2.0 m/a，从而使临邑深层地下水降落漏斗迅速扩展，成为较大的区域性降落漏斗。

4 德城区沉降中心的沉降与回弹

4.1 德城区主要沉降中心

2010年之前连续7 a的InSAR时序分析结果(图3)表明，德城区有3个典型沉降中心，分别位于德城区陈庄(德棉一厂)、长庄(华源纺织厂)和肖何庄。其中陈庄沉降中心呈南北向分布，长约2.2 km，宽1.5 km，中心最大沉降速率达到55 mm/a;长庄沉降中心呈北东向，长为2 km，宽为1.5 km，中心最大速率也超过50 mm/a。连片分布沉降出现在德城区东部的付庄和宋官屯镇，这一地区为近年来快速发展的新区。根据1991—2010年水准测量结果[11]，德城区沉降中心累计沉降量为－1 186.9～－636.9 mm，位于德棉一厂(即陈庄漏斗)位置，中心多年平均速率为59.35 mm/a。

图3 德城区2004—2010年CTInSAR监测平均地面沉降率Fig.3 Average subsidence velocity of Decheng district from 2004 to 2010 derived by CTInSAR data

其中，2005—2006年沉降量为－56 mm;2006—2007年沉降量为－89.0 mm;2007—2010年总沉降量为93.2 mm，年平均沉降量为31.1 mm，与 InSAR监测反映的德城区陈庄沉降中心的变化完全一致。

4.2 地面沉降与回弹的时间效应

由2004—2010年InSAR监测结果发现:德城区多个沉降中心年内均表现出季节性沉降与回弹特征(图4)。陈庄、长庄、肖何庄以及德棉二厂4个典型沉降中心均出现回弹，回弹最显著的位于陈庄、长庄和德棉二厂沉降中心。其中陈庄的沉降和回弹幅度均为最大，对应的年份为2007年，2—8月间沉降幅度超过80 mm，至翌年2月回弹15～20 mm。除德棉二厂外，其他3个中心2008年沉降幅度均有所减缓，范围缩小，2009年继续减缓，而市区东部在8月份之后沉降加快，2010年继续增加。

图4 德城区典型沉降中心的季节性沉降与回弹Fig.4 Seasonal subsidence and rebound of the major center in Decheng district

图4 所示5个年度的变化特征表明，以陈庄为代表的沉降中心具有明显的季节性变化特征，表现为每年的2—3月开始快速沉降，至8—9月达到最大，此后逐步回弹，至翌年2—3月重新开始沉降，沉降与回弹之比约为4∶1。

5 地面沉降及回弹与地下水位变化的关系

5.1 持续沉降与深层地下水位的关系

德州地区地下水系统分为4个孔隙含水岩组:第Ⅰ含水岩组潜水层，是主要开采层;第Ⅱ含水岩组属于咸水层，基本不开采;第Ⅲ含水岩组，底板埋深450～500 m;第Ⅳ含水岩组，底板埋深800～950 m。第Ⅲ，Ⅳ含水岩组属于承压含水层，水头埋深70～100 m，水质普遍较好，是城镇居民生活和工业用水的主要开采层，德州地下水位降落漏斗属于该层位，其中心水位(位于陈庄国棉一厂)埋深超过120 m[10]。根据德州深层地下开采引发地面沉降变化阈值的研究[11]:持续开采深层地下水是地面沉降的主要原因，地面沉降与深层水位降落呈现出显著的线性关系，深层地下水位每下降1 m，产生的地面沉降约为17.5 mm。该经验值与近4 a来地面测量结果一致，速率为30～40 mm/a。

图5所示为德城区202长观孔(位于陈庄漏斗区)1996—2010年的深层水位变化情况，表现出持续下降和年内波动的特点。德城区浅层地下水位变化曲线如图6所示。

图5 德城区202长观孔深层地下水位变化曲线Fig.5 Dynamic curve of deep groundwater level of No.202 hole in Decheng district

图6 德城区浅层地下水位变化曲线Fig.6 Dynamic curve of shallow groundwaterlevel in Decheng district

5.2 季节性沉降－回弹与地下水位波动的关系

以陈庄漏斗为例，其动态变化(图3和图4)与深、浅层地下水位的季节性波动密切相关。每年1—5月降水少，农业用水量大，浅层地下水位急剧下降;7—8月雨量增多、河道水量补给和开采量相对减少，地下水埋深逐渐回升，至年底达到次年最高点。深层承压水位变化主要取决于开采量，年水位变化与开采量关系密切。每年3—7月，深层地下水连续大量开采，漏斗内压力水头急剧下降;8月—翌年2月由于开采量减小，加之大量引入黄河水，水位有所回升，回升幅度一般小于2.0 m。

开采量决定了地面沉降的范围和强度。目前，在深层水位降落漏斗中心区(陈庄)，开采井密度为334眼/km2。2007年以前，由于深层地下水开采量大，漏斗内压力水逐年下降，年下降幅度较大，一般2～4 m;2008年以后，由于采取封井措施，深层地下水开采量减少，水位呈现回升趋势[11]。

陈庄、长庄等中心的沉降与回弹的波动与上述变化十分吻合。图7分别给出了陈庄(国棉一厂)、长庄(华源纺织厂)、国棉二厂、肖何庄和宋官屯5个沉降中心连续7 a的沉降过程(监测日期为每年的1月30日)，均表现出季节性的非线性波动(红色线表示初次估计沉降速率)。以陈庄为例，其年度内沉降与回弹与近年深层水位波动相对应，沉降速率在2008年后出现回落，小于2007年之前。

图7 德城区典型沉降中心2004—2010年地面沉降及深层水位变化序列Fig.7 Subsidence history of coherent target in the typical subsidence center of Decheng district

从地下水天然补给而言，春旱夏涝的气候特征决定了德州地区春季地下水补给量的不足[11]，而同时用水量又进一步增大，特别是陈庄、长庄和国棉二厂3个沉降中心，均为棉纺、化工等耗水企业，需要大量采水，而汛期涝季时，地表径流补给及时，地下水开采总量减小，水位出现回升，地面沉降也随之减缓，并在部分时段内出现回弹的现象。

6 结论

1)2004—2010年间的连续InSAR监测有效揭示了德州地区区域性地面沉降和典型漏斗季节性地面沉降－回弹的发展与变化特征。德州市及其周边主要地下水降落漏斗与地面沉降中心相互一致，证明了深层地下水开采是本地区地面沉降的主要原因。

2)德城区陈庄、长庄等典型漏斗的季节性沉降－回弹与该地区深浅层地下水位波动相关，地面沉降与地下水位变化表现出良好的相关关系，变化时间基本同步，滞后期约1个月。地面沉降漏斗呈现出3—8月快速沉降，沉降幅度达50～80 mm，8月—翌年2月地面回弹，回弹幅度达10～20 mm的变化特点。

志谢:感谢山东省鲁北工程地质勘察院邹祖光、陈松、杨亚宾以及河北水文地质工程地质四队杜兴明、田小伟等专家的指导和深入讨论。

[1]殷跃平，张作辰，张开军.我国地面沉降现状及防治对策研究[J].中国地质灾害与防治学报，2005，16(2):1－8.Yin Y P，Zhang Z C，Zhang K J.Land subsidence and countermeasures for its prevention in China[J].The Chinese Jounal of Geological Hazard and Control，2005，16(2):1－ 8.

[2]葛大庆，郭小方，张 玲，等.华北平原地面沉降区InSAR监测[R].北京:中国国土资源航空物探遥感中心，2010.Ge D Q，Guo X F，Zhang L，et al.Subsidence monitoring with SAR interferometry in the North China Plain[R].Beijing:China Aero Geophysical Surveying and Remote Sensing Center for Land and Resources，2010.

[3]王 艳，葛大庆，郭小方，等.长江三角洲地区地面沉降InSAR监测[R].北京:中国国土资源航空物探遥感中心，2010.Wang Y，Ge D Q，Guo X F，et al.Subsidence monitoring with SAR interferometry in the Yangtze River Delta area[R].Beijing:China Aero Geophysical Surveying and Remote Sensing Center for Land and Resources，2010.

[4]Ferretti A，Prati C，Rocca F.Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry[J].IEEE Trans Geosci Remote Sensing，2001，38(5):2202－2212.

[5]Mora O，Mallorqui J J，Broquetas A.Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images[J].IEEE Trans Geosci Remote Sensing，2003，41(10):2243－2253.

[6]Charles W，Wegmuller U，Andreas W，et al.Interferometric point target analysis with JERS－1 L－band SAR data[C]//Proceedings of the IEEE Transactfions on Geoscience and Remote Sensing.Toulouse，France，2003:4359－4361.

[7]Hoffmann J，Zebker H A，Galloway D L，et al.Seasonal subsidence and rebound in Las Vegas Valley，Nevada，observed by synethetic aperture Radar interferometry[J].Water Resources Research，2001，37(6):1551－1566.

[8]王 艳，廖明生，李德仁，等.利用长时间序列相干目标获取地面沉降场[J].地球物理学报，2007，50(2):598－604.Wang Y，Liao M S，Li D R，et al.Subsidence velocity retrieveal from long term coherent targets in Radar interferometric stacks[J].Chinese Jounal Geophys，2007，50(2):598－604.

[9]http://www.terrafirma.eu.com/documents.htm

[10]王小刚，陈 松，王秀芹，等.德州市地面沉降成因及防治对策浅析[J].地质灾害与环境保护，2006，17(3):62－66.Wang X G，Chen S，Wang X Q，et al.A Discussion on the causes of the ground subsidence and its countermeasures in Dezhou City[J].Journal of Geological Hazards and Environment Preservation，2006，17(3):62－66.

[11]杨丽芝.德州深层地下水位降落漏斗演变机制与可调控性研究[D].北京:中国地质科学院，2009.Yang L Z.Form principle and controlling－ adjusting research about deep ground water depression cone in Dezhou[D].Beijing:Chinese Academy of Geologieal Seienee，2009.

[12]王胜岭，宋 波，王德生，等.德州市临盘采油区地面沉降监测[J].山东国土资源，2009，25(1):25－28，32.Wang S L，Song B，Wang D S，et al.Land subsidence monitoring of Linpan oil mining area in Dezhou City[J].Land and Resources in Shandong Province，2009，25(1):25－28，32.

[13]何国荣，冯在敏，冯克印，等.水准测量网监测在德州市地面沉降中的应用[J].山东国土资源，2012，28(6):28－30.He G R，Feng Z M，Feng K Y，et al.Application of leveling surveying for subsidence monitoring in Dezhou City[J].Land and Resources in Shandong Province，2012，28(6):28－30.