地层形变与地下水位分层监测数据的交叉小波分析

孔祥如

(北京市地质环境监测所，北京 100195)

地面沉降是指由于自然或人为因素引发的地下松散岩层固结压缩，并导致一定区域范围内地面高程降低的地质现象，是一种缓变型地质灾害[1]，地下水超量开采是城市中发生地面沉降的主要原因[2]。L.G.Devin等[3]研究记录了中国、意大利、美国和墨西哥等近90 个国家因超采地下水造成不同程度的地面沉降。我国也有超100 座城市遭受地面沉降威胁，严重沉降区面积为7.9 万km2(累积沉降量＞200 mm)[4-5]。北京平原区自20 世纪50 年代开始出现地下水过度开采现象，导致地下水位持续下降，逐渐使北京发展成为地面沉降严重地区[6]。根据研究资料显示，北京地面沉降仍处于快速发展阶段[7]，截至2016 年，北京平原区最大累计沉降量达1 864 mm，最大年沉降速率超过110 mm/a[8]，地面沉降已经成为北京市经济社会发展的主要制约因素之一。

长期超量开采地下水导致地下水位大幅下降，邻近黏性土层向含水层释水，孔隙水压力消散，导致黏性土层压缩，形成地面沉降[9]。受压缩层的岩性、厚度和固结程度等因素控制，从水位下降到地面沉降发生这一过程具有一定滞后性[10]。对于滞后性的研究一直是地面沉降研究领域的重点。张建伟等[11]分析了上海城区地下水采灌量、开采层次及地面沉降变化特征，通过数学统计方法得出地面沉降与地下水开采量存在明显正相关，并具有时滞效应。施小清等[12]通过分析常州市分层标及水位观测孔数据并结合压缩试验，发现常州地区含水砂层在地面沉降过程中同样存在变形滞后现象，并揭示了其蠕变特性。徐进等[13]基于Biot 理论的地面沉降耦合模型，利用半解析数值原理和黏弹性流变理论，推导了可压缩土层黏弹性耦合变形的求解格式，反映土体黏滞性所导致的变形滞后效应。Zhou Zhifang 等[14]提出了释水滞后率的概念用来描述弱透水层中的滞后释水规律，其研究表明弱透水层的水动态响应滞后效应不仅与弱透水层的性质有关，也与弱透水层厚度的平方成正比关系。已有研究主要侧重于土层的滞后机理分析，对滞后时间的研究也主要是通过室内试验或数值模拟得出理论或预测结果，对实际监测数据中的准确滞后时间缺乏有效的识别方法。

小波变换是一种数据时频局部化分析工具，由其衍生的交叉小波变换能有效识别出两组序列的周期相关性，在天文、气象、地学等领域已取得一些应用成果[15-17]，在地下水学科方面，主要用于分析大气降水与地下水位的周期性变化特征[18-20]，而对于地面沉降尤其是分层沉降与地下水位的时滞响应方面缺乏相关研究应用。笔者选取北京平原区地面沉降监测站地下水位与分层标的长时间序列监测数据，采用Mann-Kendall趋势检验、连续小波变换、交叉小波变换等方法，对不同深度的地层垂向形变量和水位变化数据进行处理和对比分析，获取地层形变和地下水位动态的周期规律，继而进行交叉小波变换定量分析不同埋藏深度下的地层压缩-反弹对地下水位动态响应的滞后特征，以期进一步认识地面沉降成因机理，提高地面沉降预测模型精度。

1 研究区概况

北京平原位于华北平原北部，北依军都山，西靠太行山余脉西山。在大地构造上属于中朝准地台(Ⅰ)华北断陷拗(Ⅱ)，中生代以来，受印支、燕山、喜马拉雅多期构造活动作用，进一步形成隆起-凹陷相间的次级构造单元。北京平原第四系属典型山前冲洪积地层，整体上由永定河、潮白河、温榆河、泃河和拒马河五大水系的冲洪积扇群和多个沉积凹陷构成，第四纪沉积物岩性横向变化迅速，自山前向平原沉积物颗粒直径逐渐减小，在垂向上，则呈砾石层、砂层、黏土层交互出现的互层结构。第四系孔隙水是研究区地下水主要开采层。根据地层年代、岩性、埋藏条件及地下水补径排条件，通常将北京平原第四系含水层概化为3 组：第一组含水层(Ⅰ)为潜水和浅层承压水，对应地层为全新统(Q4)和上更新统(Q3)；第二组含水层(Ⅱ)为中深层承压水，对应地层为中更新统(Q2)；第三含水层(Ⅲ)为深层承压水，对应地层为下更新统(Q1)。每一组含水层多夹杂数层以粉质黏土、黏土或粉土为主的弱透水层(图1)。由于地下水超量开采导致含水层水位下降，弱透水层向含水层释水使得弱透水层孔隙水压力降低，进而引发黏性土层压缩，是该地区发生地面沉降的主要原因。

图1 北京平原区典型第四系地质剖面Fig.1 Typical Quaternary geological section in the Beijing Plain

2 研究方法

2.1 Mann-Kendall 趋势检验

Mann-Kendall 趋势检验是一种非参数统计检验方法，具有样本不必遵从某种特定分布、结果也不受少数异常值干扰的特点，适用于长时间序列的趋势分析，被世界气象组织推荐并广泛应用，在水文学领域也积累了大量研究成果[21-24]。

设xn(n=1,2,···,N，N为样本数量)为平稳独立序列，统计量U定义为：

当n≥10，统计量U近似服从正态分布，不考虑序列中等值数据点情况，其均值E(U)=0。

标准化的检验统计量ZMK采用下式计算：

式中：g、h为任意样本序号，且g＞h。当β为正值，表示上升趋势；当β为负值，表示下降趋势[19]。

2.2 连续小波变换

小波变换将波频显示在不同时间尺度上，来反映长时间序列中不同尺度显著周期及其所在的时间段。本文选用在时间和频率的局部化上具有良好平衡性的Morlet 小波作为母小波φ0，函数为：

式中：ω0为 无量纲频率；η为无量纲时间。

时间序列的连续小波变换可定义为对序列xn基于小波标准化算法的卷积，函数为：

式中：δt为均一时间步长；s为小波尺度。

背景功率谱采用红噪声检验，红噪声检验过程用一阶自回归方程。背景红噪声功率谱Pk定义为：

式中：r为红噪声功率谱中自回归方程的相关系数；k为傅里叶频率系数。

2.3 交叉小波变换

交叉小波变换能够分析两组序列在时频空间中能量共振和协方差分布规律，可以在不同时间尺度上揭示两组序列频率周期的一致性和相关性，并能通过小波相位角分析，计算时频空间中的相位关系。

式中：σX、σY分别为时间序列Xn和Yn的标准差；当背景功率谱为实小波时，v=1，为复小波时v=2；Zv(p)是概率p的置信度水平。

交叉小波的相位角am定义为：

式中：ai(i=1,2,···,N)为小波变换时间域上的单个相位角。

小波相位角的标准偏差为：

本文所采用的小波分析计算方法和程序代码主要参考C.Torrence 等[25]和A.Grinsted 等[26]的研究成果。

3 数据来源与预处理

3.1 数据来源

本文使用的地面沉降和地下水水位数据来自北京平各庄地面沉降监测站，该站建于2008 年，站内埋设有分层标和地下水位监测井，监测层位包含北京平原区第四系主要含水层(砂层)和弱透水层(黏性土层)，地层结构具有代表性。本文使用的数据时间段为2008 年6 月-2018 年12 月，监测频率为1 次/d。地下水位监测井(D-1-D-5)主要监测不同含水层水位，地下水位时间序列(D1-D5)为单井监测数据；分层标监测井(F-1-F-5)监测2 个分层标之间的地层，包括弱透水层和含水层，其地层形变时间序列(F1-F5)为地层顶底2 个分层标监测数据之差，即监测层位的垂向形变量。地下水位和沉降监测层的基本情况见表1。

表1 地下水位和沉降监测层基本情况Table 1 Basic information of groundwater monitoring wells and subsidence monitoring layers

3.2 数据预处理

地下水长期持续开采、地应力载荷等长时间尺度因素会对数据分析产生干扰，因此首先对地下水位和地层形变数据进行Mann-Kendall 趋势检验，以判定其变化趋势。结果见表2，可以看出，各时间序列的ZMK值均为负值且均小于-3.3，为下降趋势，ZMK通过了α=0.001水平的显著性检验，下降趋势显著。地下水位方面，D-1 和D-2 井监测层位的深部承压含水层为地下水主要开采层，2008-2018 年地下水位均呈持续下降趋势，地下水位年均下降速率分别为1.780 3 m/a和1.184 9 m/a；D-3 中部承压含水层在2014 年以前水位持续下降，2014 年以后受地下水禁限采政策影响开采量逐渐减少，该含水层水位出现缓慢回升，2008-2018 年地下水位年均下降速率为0.252 1 m/a；D-4 和D-5 井监测层位的浅部含水层地下水开采量较小，主要受大气降水影响，10 a 平均下降速率分别为0.416 7 m/a和0.491 7 m/a。地层形变方面，F-1 地层自监测以来持续发生压缩沉降，由于下更新统土体固结程度较高，可压缩性较低，属于地面沉降次要贡献层，10 a 平均沉降速率为1.772 8 mm/a；F-2 与F-3 地层是地面沉降的主要贡献层，两者沉降量占地面总沉降量的80%以上，其中F-2 地层10 a 平均沉降速率为13.046 9 mm/a，2014 年以后受地下水位回升影响，F-3 地层沉降速率有所减缓，但仍然保持沉降，10 a 平均沉降速率为

表2 Mann-Kendall 趋势检验 ZMK 和 β值Table 2 ZMK and β of the Mann-Kendall trend test

7.158 7 mm/a；F-4 和F-5 地层属于弱沉降层，地层沉降量较小，个别年份还会有所反弹，10 a 平均沉降速率分别为0.866 4 mm/a 和0.709 0 mm/a。为了消除变化趋势的影响，根据R.T.Hanson 等[27]的数据前处理方法，将地下水位和地层形变数据减去三次多项式回归方程，取其残差作为分析地下水位和地层形变关系的时间序列数据。

4 结果与分析

4.1 连续小波变换

将全部的地下水位和地层形变量时间序列分别进行连续小波变换(图2)，以分析其周期性波动特征和显著周期时段等信息。

地下水位方面(图2a)，埋深81 m 以下的中-深层承压水时间序列D1、D2 和D3 存在1 a 左右的主震荡周期，连续分布且周期显著；60 m 以浅含水层地下水位在大部分时域上周期性不显著，其中浅层承压水时间序列D4 仅在2008-2012 年存在1 a 左右的震荡周期，潜水时间序列D5 仅在2014-2015 年存在0.15～0.75 a 的主震荡周期。2000 年以后，由于浅部含水层水位下降和水质恶化，北京平原区地下水开采深度不断增加，中-深层承压水成为主要开采地下水[28]，尤其在农作物灌溉期地下水超采现象严重，造成深部含水层水位具有明显的季节性波动特征[29]，是深部含水层较浅部含水层水位震荡周期更为显著的主要原因。

图2 各时间序列连续小波变换Fig.2 Continuous wavelet transforms of the time series

地层形变方面(图2b)，根据监测数据埋深在63～209 m 的地层沉降占比在80% 以上，属于沉降主要贡献层，其中地层形变时间序列F2 和F3 在全时域存在1 a 左右的主震荡周期，且通过95%红噪声检验的区域与同层水位具有相似形态，表明该严重沉降层的地下水位波动可能与地层形变量存在时频域相关；埋深超过209 m 的地层形变量相对较小，其时间序列F1 在2010-2017 年存在1 a 左右的主震荡周期，通过95%红噪声检验的区域与同层水位较为相似；埋深小于63 m 的地层形变量最小属于弱沉降层，其中时间序列F4 在2010-2017 年存在1 a 左右的主震荡周期，时间序列F5 在2013-2015 年存在0.5～3.0 a 的主震荡周期，与同层水位小波变换形态相比均具有更显著的周期性。

4.2 交叉小波变换

对相同层位的地下水位和地层形变量时间序列进行交叉小波变换(图3)，以分析其共振周期、显著时段和相位关系等，研究地面沉降对地下水位变化的时滞响应特征。

由图3a-图3c 可知，埋深在63 m 以下的地层地下水位和地层形变量通过95%红噪声检验且共振周期在0.75～1.29 a 的区域能量密度最高，显著时段连续分布整个序列时域，表明两者共振关系显著。滞后性方面，埋深63～209 m 地层形变量相对地下水位的迟滞时间相近，分别为(7.16±7.09) d 和(9.66±6.62) d(表3)；埋深209 m 以下地层形变滞后时间更长，为(16.58±8.91) d。由此可以看出，深部主要沉降层的地层形变与地下水位动态相干性较强，形变周期相对地下水位周期具明显的滞后性，但不同层位滞后时间具有差异性。

图3 地下水位与地层形变量交叉小波变换Fig.3 Cross wavelet transforms of groundwater levels and stratum deformation values

由图3d 和图3e 可知，埋深在63 m 以浅的地层地下水位和地层形变量通过95%红噪声检验的区域在时频域上不连续，在部分时域上具有显著共振周期。埋深32～63 m 地层形变滞后时间为(21.01±21.06) d(表3)，相较于深部主要沉降层有更长的滞后时间；埋深2～32 m 地层交叉相位的相位误差值远大于平均相位值，即在共振周期显著的时频域里形变滞后时间普遍存在负值，说明浅层地层的形变量与潜水水位变化不具有可信的时滞关系。

表3 地下水位和地层形变量交叉小波变换周期与相位特征Table 3 Periodicities and phases of cross-wavelet transforms of groundwater levels and stratum deformation values

为进一步验证埋深63 m 以浅地层形变是否与不同层位含水层地下水位存在迟滞关系，对埋深63 m 以浅地层的地下水位和形变量进行错层交叉小波变换(图4)。

图4 浅部错层地下水位与形变量交叉小波变换Fig.4 Cross wavelet transforms of groundwater levels and stratum deformation values in shallow staggered layers

图4a 看出，埋深在81～119 m 含水层地下水位与埋深在32～63 m 地层形变量在全时域具有0.74～1.25 a的显著共振周期，形变时滞为(32.02±9.67) d(表4)，说明中层承压水的水位波动对上部邻近地层的形变也具有一定影响，但形变滞后时间较长，引发的形变量较小；图4b-图4c 看出，埋深在60 m 以浅地层的含水层地下水位与相邻地层形变量通过95%红噪声检验的区域形变时滞分别为(3.90±64.87) d 和(40.87±57.38) d(表4)，表明潜水和浅层承压水与相邻地层的形变不具有明确的时滞关系。

表4 浅部错层地下水位和形变量交叉小波变换周期与相位特征Table 4 Periodicities and phases of cross-wavelet transforms of groundwater levels and stratum deformation values in shallow staggered layers

5 结 论

a.研究区作为地下水主采层的中-深层承压水含水层水位具有1 a 左右的主震荡周期，在全时域连续分布且周期显著，符合地下水位季节性波动的规律特征；潜水和浅部承压水含水层水位仅部分年份通过95%红噪声检验，在大部分时域无显著周期。

b.深部严重沉降层的形变量具有1 a 左右的主震荡周期，在全时域连续分布且周期显著，地下水位与形变量通过95%红噪声检验的共振周期为0.75～1.29 a且连续分布，验证了两者具有可靠相关性，地层形变滞后于水位变化，地层由浅到深形变滞后时间分别为(16.58±8.91) d、(7.16±7.09) d 和(9.66±6.62) d。

c.浅部弱沉降层中，埋深在2～32 m 地层形变量在COI 区域内存在0.5～3.0 a 的主震荡周期，与地下水位不具有时滞关系；埋深在32～63 m 地层形变量具有1 a 左右的主震荡周期，与第一组含水层水位不具有时滞关系，而与中层承压水含水层水位存在显著共振周期，形变时滞为(32.02±9.67) d。

d.运用小波变换方法分析北京平各庄地面沉降监测站不同深度地层的形变滞后特征，精细刻画地面沉降与地下水的相互关系，为构建地面沉降精细化模型、提高地面沉降预测精度以及研究更有效的地面沉降防控措施提供新的技术思路。但地面沉降的发展受岩性、厚度、地应力和固结程度等多种因素的控制，地面沉降滞后性研究还需大量数据的支撑。后续将开展其他监测站的数据分析，扩大时间尺度，用更充足的数据去论证地面沉降对地下水的响应特征；另一方面将加强对研究区钻探、物化探及相关试验数据的搜集，结合小波变换分析结果，进一步研究不同环境因素对地面沉降滞后性的影响。