0 引言

地面沉降是在自然和人为双重因素的影响下，由于地壳表层土体的压缩而导致区域性地面标高缓慢降低的一种不可补偿的永久性环境地质灾害。地面沉降具有生成缓慢、持续时间长、影响范围广、成因机制复杂和防治难度大等特点^[1]。从20世纪60年代以来，北京平原区地面沉降一直呈快速发展的趋势。早在2003年，北京平原出现了昌平沙河—八仙庄、顺义平各庄、东北郊—来广营、东郊八里庄—大郊亭和大兴榆垡—礼贤五大沉降区，形成了近1 000 km²的地下水位降落漏斗。研究^[2]表明，北京市地面沉降的主要原因是长期的地下水过量开采。近年来，随着南水北调工程以及地下水限制开采政策的实施，北京市地下水开采量有所减缓。但据统计，目前北京市50%~60%的供水仍来自地下水，由此引发的地面沉降仍未得到有效控制。因此，通过地面沉降监测及其周期性分析，掌握地面沉降演化规律，对制定控制地面沉降、限制地下水开采措施有重要的支撑作用。

传统的地表形变沉降监测手段主要包括水准测量和GPS。然而水准测量方法耗时长，监测范围有限；GPS测量设备造价昂贵，定期维护费用高，监测密度小。相比之下，合成孔径雷达干涉(InSAR)技术可以获取大范围、高精度的地表形变信息^[3]。为了克服时间基线过长造成的失相干，Ferretti等^[4]在1999年提出了永久散射体干涉测量(persistent scatterer Interferometry, PSI)技术。PSI技术选择高相干点的永久散射体，极大地提高了长时间序列干涉测量的精度和可靠性。干涉点目标分析(IPTA)方法由Charles Werner^[5]于2003年提出。其核心思想是：对具有稳定频谱特性的点目标相位信息进行时间维和空间维的复合分析，得到其时序形变信息。IPTA技术同样可以克服时间基线过长引起的失相干，同时提高数据处理效率^[6-8]。

目前，国内外学者针对地面沉降的演化特征、成因机理等方面有许多讨论和研究。在演化特征方面，研究者^[9-10]结合水文地质条件、土地利用类型差异等研究地面沉降的演化规律。葛大庆等^[11]利用InSAR长时间序列研究德州市地面沉降-回弹与地下水波动特征，结果表明地面沉降与地下水位变化密切相关。雷坤超等^[12]基于PS-InSAR技术建立地面沉降综合分析模型，揭示了北京区域地面沉降的不均匀性，发现北京的多个沉降漏斗连成一片，沉降受雨水渗入补给影响较大，有季节波动性特征。成因机理方面的研究^[13-14]表明，长期超量开采地下水和天然气等资源是目前平原区地面沉降的主要诱发原因。例如，贾三满等^[2]认为北京平原区地面沉降主要是由于超量抽汲地下水引起。Estelle Chaussard等^[15]发现，地下水及天然气的开采是造成印度尼西亚地面沉降的主要原因。陈蓓蓓等^[16]采用InSAR技术、多源遥感技术与水文地质学交叉研究，揭示了北京平原区地下水演化与地面沉降响应机理；还发现北京平原区地面沉降与建筑载荷密度存在一定的相关关系。付延玲等^[17]研究发现由高层建筑荷载引起的地面沉降是以建筑物中心为漏斗中心的漏斗状。周超凡等^[18]采用数据场模型发现了北京市地面沉降较为严重的朝阳区北部与中部区域有高交通载荷。

小波分析是一种信号时、频局部化分析的新方法，在信号处理、图像压缩等领域有广泛的应用^[19]。其原理是通过将时间序列分解到时间频率域内，可以得到时间序列的显著波动模式，即周期变化动态。Morlet小波变换能有效地分离出隐藏在时间序列中的多种变化周期。早在2010年，研究者^[20-22]将小波分析引入地面沉降研究，拓宽了其应用范围。最先将小波分析引入地理学的是Grinsted等^[23]，其使用交叉小波变换(XWT)和小波相干性(WTC)研究北极涛动与海冰范围特征。在国内，朱锋等^[24]将小波方法应用到了地面沉降研究中，用bd5小波识别京津冀高铁北京段的不均匀沉降段；高明亮等^[25]针对北京国际机场地面沉降问题，采用了连续小波变换(CWT)分析地面沉降与地下水的关系。

本文选取27景Radarsat-2数据，采用IPTA方法获取北京平原区域地面沉降监测信息，结合Morlet小波分析方法，分析了监测井的地面沉降周期性演化特征。

1 研究区概况

研究区位于北京市朝阳区东部、顺义区南部和通州区西部的交界区域，地理位置为116°24′00″E—116°45′00″E，39°48′00″N—40°07′00″N，面积为731.50 km²(图 1)。北京市由西北向东南、从山前至平原，第四系沉积厚度逐渐增大、层次增多，沉积物颗粒变细。研究区第四系黏性土可压缩层厚度为60~250 m^[26]。自1987年至今，北京市地下水持续超量开采，使北京平原区的地下水位逐年下降，导致土层压缩，从而引发了地面沉降。目前，区域内包括东郊八里庄—大郊亭、东北郊—来广营2个地面沉降区；且地面沉降呈逐年加剧，沉降漏斗连成一片的趋势。本文在研究区中任意选择4个监测井，分别为1#、2#、3#、4# (图 1)，将距离各监测井200 m内永久散射体(persistent scatterer，PS点)沉降均值作为该点沉降值，研究其时间序列周期性沉降特征。

2 研究方法 2.1 基于IPTA的监测方法

IPTA方法是Charles Werner^[5]于2003年提出的，其利用迭代回归分析来获取形变历史、高程改正值、大气效应值等，具体步骤包括：

1) PS点的提取。选择主影像，配准合成孔径雷达(synthetic aperture radar，SAR)影像集，裁剪研究区，根据光谱属性或者稳定的后向散射强度选择PS候选点。

2) 生成初始差分干涉图。根据主影像生成点差分干涉图；对于每一个干涉像对，用外部数字高程模型(digital elevation model，DEM)模拟研究区PS点高程相位，并从干涉相位中减去模拟的高程相位，得到差分干涉相位φ_dif。

3) 对PS点差分干涉相位进行回归分析获取线性形变信息。在IPTA的相位模型中，φ_dif被分解为地形误差相位φ_topoerr、形变相位φ_def、大气效应相位φ_atm和噪声相位φ_noise，即

回归分析过程中将得到线性形变速率、高程改正值、残余相位。其中形变相位φ_def分为线性和非线性2个部分，非线性形变暂时归入残余相位。

4) 非线性形变提取。回归分析得到的残余相位包括大气效应相位、非线性形变相位、基线误差和噪声相位等；之后，根据不同相位成分在时间和空间上表现出的相关性差异，分别进行时间域和空间域的滤波处理，从而分离大气相位和噪声，得到非线性形变相位。

5) 提取形变信息。将上一步骤中得到非线性形变融合到线性形变速率中，并重新计算相干点形变速率。通过对各时段沉降速率在时间上进行积分处理，获得对应每个时间点的累积形变量。

考虑到研究区和现有数据覆盖情况，本次研究选择的是2011—2014年的27景RadarSat-2数据。RadarSat-2是加拿大太空署与加拿大麦克唐纳·迪特维利联合公司(MacDonald Dettwiler and Associates Ltd.，MDA)合作的商用雷达数据。本研究所用到的数据产品，其距离向分辨率为11.8 m，方位向分辨率为5.2 m，波长为5.6 cm(为C波段数据)，重访周期为24 d。用于去除地形相位和地理编码的外部数据为美国航天飞机测图任务(SRTM3)获取的90 m分辨率的DEM。

2.2 Morlet小波分析方法

20世纪80年代初，由Morlet提出的一种具有时-频多分辨功能的小波分析(wavelet analysis)为更好地研究时间序列问题提供了可能。它能清晰地揭示出隐藏在时间序列中的多种变化周期，反映系统在不同时间尺度中的变化趋势，并能对系统未来发展趋势进行定性估计^[19]。本次研究采用Morlet小波为基小波，运用连续小波变换分析北京平原地面沉降周期性特征。

对于f(t)∈L²(R)，连续小波变换(CWT)被定义为信号f(t)与小波基函数的内积:

式中：f(t)为原始信号；t为时间；L²(R)为f(t)的矢量空间；W_f(a, b)为小波变换函数；a为尺度因子；b为定位因子；ψ_{a, b}为基小波；是复共轭函数。

将由小波变化方程计算得到的小波系数的平方值在时间域上积分，其结果即为小波方差：

小波方差随尺度a的变化过程，称为小波方差图。由式(3)可知，它能反映信号波动的能量随尺度a的分布。因此，小波方差图可用来确定信号中不同尺度扰动的相对强度及其对应的时间尺度，即主周期。

小波变换的数据要满足正态分布，运用K-S检验可知，1#、2#、3#、4#正态分布检验结果分别为0.996、0.402、0.966和0.992。其P值(K-S检验的显著性差异)均大于0.05，通过K-S检验(置信度设为95%)，沉降值时间序列符合正态分布，可以进行小波变换。

3 结果与分析 3.1 IPTA结果及精度验证

利用IPTA获取地表形变信息，在研究区范围内，共识别191 706个PS候选点。在2011—2014年3年间，研究区各PS点年均沉降速率分布如图 2a所示。从图 2a可以看出，研究区最大沉降速率为162.70 mm/a。研究区内地面沉降分布不均，该时段内累积沉降量最大值为645.86 mm。采用反距离加权(IDW)方法对PS候选点的形变速率进行空间插值，得到研究区的地面沉降速率插值结果(图 2b)，获取研究区沉降速率空间演化趋势。

采用自然间断点分类方法，将PS候选点分为7类，统计每一类的面积比例和PS点数。由表 1可知：第1类地面沉降速率(0.00~ 22.79 mm/a)及第2类地面沉降速率(22.79~37.81 mm/a)占比较大，分别为为22.44%、24.40%，包含PS点数分别为52 809、42 917；第7类地面沉降速率(111.52~162.70 mm/a)所占面积比例最小，为5.08%。总之，地面沉降速率的分布不均匀性较大，即不均匀沉降演化趋势明显。

以水准点地理位置为参考，距离水准点200 m范围内的PS点均值作为该点形变值，计算水准观测与IPTA反演结果的误差(表 2)，检验IPTA技术的监测精度。为了与IPTA获取的年均沉降速率进行比较，分别求取7个水准点在2011、2012年的年累积沉降量。校验结果显示，误差最大值为9.586 mm/a，最小值为0.906 mm/a，形变误差在1 cm/a内。在置信度为95%的条件下，相关系数达到了0.972 8，说明IPTA反演结果与水准观测结果一致。

3.2 基于Morlet小波分析典型区地面沉降周期性

小波分析要求输入的信号必须是等时间间隔的。因此，必须对现有的IPTA时间序列沉降结果进行时间域的插值。首先以4个监测井为圆心，做距离为200 m的缓冲区，缓冲区内PS点沉降均值作为该点沉降值；之后对地面沉降时间序列进行3次样条插值和误差计算，最终得到2011年6月26日到2014年11月1日期间，每隔24 d间隔的地面沉降时间序列。

将等时间间隔的沉降时间序列信号输入Morlet小波变换工具进行处理，得到小波系数，同时计算小波系数的实部、模、模方、方差，之后将小波系数方差绘制成图(图 3)。如前文所述，小波系数方差图可以清楚地展示出信号中不同尺度扰动的相对强度及其对应的时间尺度，即主周期。图 3a中，1#点小波方差图呈现出3个峰值，时间尺度分别对应为28、16、7 T(1 T表示1个24 d的时间段)。其中：28 T的时间尺度对应的峰值最大，说明28 T时间尺度处周期震荡最强，为该位置沉降的第一主周期；16 T时间尺度峰值次之，为1#点地面沉降时间序列的第二主周期；7 T时间尺度对应最小峰值，为该点地面沉降时序演化的第三主周期。上述3个周期的波动控制着该区域地面沉降在整个时间域内的变化特征。相似地，如图 3b所示，2#点存在4个峰值，时间尺度分别为28、22、13、7 T。其中22 T对应该位置的第一主周期；28、13和7 T处分别对应其第二、第三和第四主周期，4个周期的波动控制着该区域地面沉降在整个时间域内的变化特征。图 3c中，3#位置在28 T处存在1个明显峰值，其为该点地面沉降时间序列的主周期。图 3d中：4#位置在28 T处存在峰值，其为对应该点地面沉降时间序列主周期；7 T为第二主周期，2个周期的波动控制着该区域地面沉降在整个时间域内的变化特征。

根据小波方差结果，分别绘制1#、2#、3#、4#位置地面沉降演化的各时间尺度Morlet小波变换实部过程线(图 4-7)。

如图 4所示，1#在28 T的时间尺度上，地面沉降演化经过了约3个整周期，平均周期为13.3月；在16 T的时间尺度上，平均周期为7.6月；在7 T的时间尺度上平均周期为3.8月。类似地，如图 5所示：2#点在22 T的时间尺度上，地面沉降演化经过了约3.75个周期，平均周期为10.7月；在28 T的时间尺度上，平均周期为13.3月；在13 T的时间尺度上平均周期为6.7月；在7 T的时间尺度上平均周期为3.5月。如图 6所示，3#点在28 T的时间尺度上经过了约3个周期，平均周期为13.3月。如图 7所示：4#位置在28 T的时间尺度上，地面沉降演化经过了约3个周期，平均周期为13.3月；在7 T的时间尺度上平均周期为3.8月。

4 结论与建议

1) 2011—2014年，北京地区地面沉降发展迅速，研究区最大年沉降速率达到162.70 mm/a，最大累积沉降量为645.86 mm。

2) 以研究区内1#—4#为对象，通过Morlet小波变换分析不同时间尺度下地面沉降的周期性特征。结果表明，地面沉降时序演化在空间分布上呈现差异性特征。地面沉降周期随着研究尺度的不同而发生相应的变化，这种变化一般表现为小时间尺度的变化周期往往嵌套在大尺度的变化周期之中。地面沉降在时间域中存在多层次的时间尺度结构和局部变化特征。

3) 长期超量开采地下水是北京地面沉降的主要原因。作为地下水的主要补给源，降雨量的年际变化可能影响地下水位的变化，从而也可能诱发地面沉降的季节性波动。这次研究揭示了区域地面沉降多尺度时序演化特征，在未来的研究中，需要进一步研究地面沉降与地下水时序演化的关联分析。