北京南部沉降区黏土层变形特征和参数分析.pdf

高 校地质学报Geological Journal of China Universities2023 年 8月，第 29 卷，第 4 期，590-599页August 2023，Vol.29，No.4，pp.590-599 收稿日期：2021-07-30；修回日期：2021-08-22基金项目：本文为国家自然科学基金项目（批准号：41877180）；北京市财政项目（PXM2019_158305_000012）联合资助作者简介：尤斌昊，男，1995年生，硕士研究生，主要从事地下水和地面沉降数值模拟的研究；E-mail: *通讯作者：叶淑君，教授；E-mail:北京南部沉降区黏土层变形特征和参数分析尤斌昊1，叶淑君1*，田 芳2，罗 勇21.南京大学 地球科学与工程学院，南京 210023；2.北京市水文地质工程地质大队，北京 100195 摘要：北京地面沉降分南北两个大区，南部沉降区主要为礼贤榆垡沉降区，近年该区沉降速率加快，20172019年沉降速率达5075 mm/a，沉降呈现新的特征。文章基于该沉降区唯一的沉降观测站（榆垡沉降监测站）20112017年7年长时间序列的沉降和水位观测资料，首先系统分析了该地区地面沉降监测孔监测地层尤其是黏土层的变形特征，然后运用应力应变图解法和配线法计算该区黏土层的弹性贮水率、非弹性贮水率、垂向渗透系数和变形滞后时间。结果表明，各土层既有弹性变形又存在塑性变形，变形均以塑性变形为主，有明显的变形滞后性；但各土层在变形速率上呈现不同特征。图解法结果表明，第一和第二黏土层非弹性贮水率在1.4910-52.1010-5之间，弹性贮水率在3.9310-58.3810-5之间，非弹性贮水率为弹性贮水率的24倍，垂向渗透系数在6.2010-64.3910-5 m/d之间，土层变形滞后时间为2.112.8年。配线法估算出研究区浅部土层的贮水率在4.810-4左右，深部土层的贮水率在1.0810-42.0510-4之间，浅部土层的贮水率是深部土层的24倍，浅部土层垂向渗透系数比深部土层大12个数量级。关键词：地面沉降；黏土层；变形特征；水文地质参数；图解法；配线法中图分类号：P641 文献标识码：A 文章编号：1006-7493（2023）04-590-10Analysis of Deformation Characteristics and Parameters of Clay Layers in South Subsidence Area of Beijing YOU Binhao1，YE Shujun1*，TIAN Fang2，LUO Yong21.School of Earth Sciences and Engineering,Nanjing University,Nanjing 210023,China;2.Beijing Institute of Hydrogeology and Engineering Geology,Beijing 100195,ChinaAbstract:The land subsidence in Beijing is divided into north and south areas,and the southern area mainly includes Lixian-Yufa subsidence zone.The subsidence rate in this area has accelerated in recent years and reached 50-75 mm/a during 2017-2019,showing new subsidence characteristics.Based on the subsidence and water level observation data of 7-year time series from 2011 to 2017 at Yufa Subsidence Monitoring Station,the only subsidence observation station in the south subsidence area.The deformation characteristics of soil layers,especially clay layers,are systematically analyzed.Then the elastic and inelastic specific storage,vertical hydraulic conductivities and deformation lag time of the clay layers are estimated by the stress-strain diagram method and the type curve method.The results show that there are both elastic deformation and plastic deformation in each soil layer,and the plastic deformation is the main deformation with obvious deformation hysteresis.However,the deformation rate of each soil layer presents different characteristics.The results of the stress-strain diagram method show that the inelastic specific storage of the first and second clay layers is between 1.4910-5 and 2.1010-5,the elastic specific storage is between 3.9310-5 and DOI:10.16108/j.issn1006-7493.2021083引用格式：尤斌昊，叶淑君，田芳，罗勇.2023.北京南部沉降区黏土层变形特征和参数分析J.高校地质学报，29（4）：590-599尤斌昊591尤斌昊等：北京南部沉降区黏土层变形特征和参数分析4 期8.3810-5,the inelastic specific storage is 2-4 times of the elastic specific storage,the vertical hydraulic conductivity is between 6.2010-64.3910-5 m/d,and the soil deformation lag time is 2.1-12.8 years.The specific storage of the shallow layers in the study area estimated by the type curve method is about 4.810-4,and that of the deep layer is between 1.0810-4 and 2.0510-4.The specific storage of shallow layer is 2-4 times of that of deep layer,and the vertical hydraulic conductivities are 1-2 orders of magnitude larger than that of the deep layers.Key words:land subsidence;clay layer;deformation characteristics;hydrogeological parameter;stress-strain diagram method;type curve methodCorresponding author:YE Shujun,Professor;E-mail: 1 引言过量开采地下水引起的地面沉降灾害在世界很多国家和地区发生。据联合国教科文组织地面沉降国际倡议 19 位研究人员 2021 年发表在Science期刊的最新研究表明；尽管在过去的一个世纪里，有包括中国在内的 34 个国家/地区报道了过量开采地下水引起的地面沉降区 200 多处，但实际发生地面沉降的区域远大于此；预测 2040 年全球潜在沉降面积增加 7%，影响 16 亿居民（Herrera-Garcia et al.,2021）。华北平原是目前中国地面沉降范围最广、沉降速率最快、危害最为严重的地区（Guo,2015;Ye et al.,2015）。截止 2015 年，华北平原沉降区面积达 8.57 万 km2，其中以京津冀地区表现最为突出（何庆成等，2006a,b;王云龙等,2018;李文鹏等，2021）。北京至 2016 年，最大累计沉降量达 1864 mm，最大年沉降速率超过 110 mm/a，地面沉降发育程度中到强区（累计沉降量大于 300 mm 或近 5 年沉降速率大于 10 mm/a）面积达 3405 km2，约占北京平原区面积的 55%（叶超等，2019）。北京平原区地面沉降发育南北分区明显，南部沉降区主要为大兴榆垡礼贤沉降区（周毅等，2016）（图 1），该沉降区形成于 20 世纪 70 年代后期，19552012 年最大累计沉降量达1137 mm，最大年沉降速率达 78.7 mm/a，目前仍处于快速发展时期（何庆成等，2006a,b;郭海朋等,2021），20172019 年沉降速率达 5075 mm/a（晏霞等，2021），表现出沉降速率加快，深层沉降贡献增加等特征（周毅等，2016;曹群等，2019;田芳等,2017;杨艳等，2021）。北京平原区地面沉降发展与黏性土可压缩层厚度显著相关（雷坤超等，2016;卢毅等，2019;赵亚丽等，2021）。同一地区，黏性土层单位沉降量是砂层、砂砾石层的 23 倍，地面沉降主要发生在黏性土层中（叶超等，2018）。大兴榆垡-礼贤沉降区地下水超采严重，开采量超过年平均可开采资源量 20%50%（叶超等，2018）。包含地下水流模型和沉降模型的区域地面沉降数值模型是实现地面沉降科学管控的手段（叶淑君等，2005;薛禹群等，2006;李文运等，2012;熊小锋等，2017;罗跃等，2018;Ye et al.,2018;杨蕴等，2019）。北京先后建立了多个地面沉降模型（崔亚莉等，2003;杨勇等，2010;孟世豪等，2021）。以往所建立的模型采用的弹性贮水率、非弹性贮水率、垂向渗透系数等参数取值多以经验值为准，导致模型存在一定的误差。土层变形特征是地面沉降机理研究的重要组成部分，也是构建地面沉降数学模型的基础。构建研究区地面沉降数值模型需要先厘清该区域黏土层变形特征（刘妍君等，2021），进而获取地面沉降模型涉及的水力学和土力学参数（董成志等，2017;孟世豪等，2021）。上海、苏锡常和天津地区地面沉降数据监测开始时间较早，前人根据这些地区的水位和沉降监测资料对土层变形特征和力学参数做了大量研究（叶淑君等，2005;张云等，2006;施小清等，2006;于军等，2007;罗跃等，2015;郭海朋等，2017）。已有学者对北京中东部地区的土层变形特征进行分析研究（周毅等，2016;杨艳等，2021;雷坤超等，2016;田芳等，2012;罗勇等，2011），但对北京南部沉降区的黏土层变形特征，尤其是水文地质参数研究甚少。构建地面沉降模型时最大的工作量是识别模型参数，而一个好的初始参数值可以很大程度降低模型参数识别的工作量。含水层的初始参数往往可以通过野外抽水试验获得较好的初始值，但抽水试验不适用于黏性土层。而黏性土层的水文地质参数对于构建反映高校地质学报2 9 卷 4 期592实际沉降特征的模型至关重要（罗跃等，2018）。黏性土层是典型的弱透水层，获取弱透水层水文地质参数的方式除了室内实验和模型反演（van der Kamp et al.,2001;Yu et al.,2013），还有根据长时间序列的水位和沉降观测资料进行计算的 图 解 法（Riley,1969;Burbey,2001）和 配 线 法（Wikramaratna et al.,1985）。室内实验方法受尺度效应影响，模型反演方法所需资料多且反演过程计算量大。应力-应变图解法可以根据沉降和相应水位资料计算土层弹性贮水率、非弹性贮水率，垂向渗透系数和变形滞后时间。特别是利用野外长时间序列水位和沉降观测资料获取的水文地质参数更能代表野外尺度含水层水力学和土力学特性，因而也更适用于区域地面沉降模型。叶淑君等（2005）应用图解法获取了上海地区弱透水层水文地质参数。庄超等（2017）提出的配线法可计算弱透水层的垂向渗透系数和贮水率，应用该方法确定了上海F10-7 分层标处第 2 弱透水层的渗透系数，此方法能够反映弱透水层释水变形过程的滞后性。北京南部沉降区唯一的沉降观测站榆垡沉降监测站，拥有 20112017 年 7 年长时间序列的沉降和水位观测资料，为应用图解法和配线法获取该地区弱透水层水文地质参数提供了宝贵数据。因此本文根据北京榆垡监测站 7 年（2011 年1 月至 2017 年 12 月）的长时间序列地面沉降和水位观测资料，在分析南部沉降区地下水引起的地面沉降过程各黏土层的变形特征基础上，计算北京南部沉降区各黏土层水文地质参数，进一步揭示南部区沉降特征和参数特性，为建立该地区地面沉降数值模型，为确定地面沉降防控措施提供参考。2 不同深度土层变形特征2.1 地层结构北京南部沉降区属永定河地下水系统，由于河图1 北京市平原区19552012年累计地面沉降量图（据周毅等，2016重绘）Fig.1 Cumulative land subsidence in plain area of Beijing from 1955 to 2012593尤斌昊等：北京南部沉降区黏土层变形特征和参数分析4 期流冲积、洪积作用的不同，堆积物颗粒组成、胶结程度、结构构造存在很大差异，从而造成岩性特征的巨大差别。大兴西北部为永定河近期多期古河道的主体部位，第四纪地层发育齐全，厚度由西南向东北逐渐增大，西部约 40 m，东部则大于 100 m。大兴区中部为永定河多期古河道扇中所在，地层岩性复杂。大兴区东南部第四纪地层发育齐全，和下伏新近系为整合接触，沉积厚度由北向南、自东向西逐渐增大（张安京等，2008;叶超等，2018）。大兴榆垡监测站监测 334 m 以浅第四系松散沉积层的水位和变形。根据区域水文地质资料，概化为 5 个含水层组，4 个黏性土层；在不同深度设 7个分层标和 5 个水位监测井。7 个分层标监测深度分别为 2 m、27 m、53 m、82 m、116 m、170 m 和310 m（图2），5个水位监测井监测深度分别为34 m、83 m、117 m、178 m 和 310 m（刘予等，2007）。其中黏性土层厚度总计 109.25 m，占地层总厚度 32.7%，2011 年 1 月至 2017 年 12 月黏性土层累计沉降 120.64 mm，占总沉降量的 49.5%，单位厚度黏性土层累计沉降量（1.81 mm/m）约为含水层（0.91 mm/m）的 2 倍。不同深度含水层 2011 年 1 月至 2017 年 12 月水位和累计变形量随时间变化的关系如图 3图 5所示，其中黏性土层水位根据上下含水层水位插值所得（庄超等，2017），水位变化对应于黏性土层中有效应力变化，累计沉降量对应于土层应变，作应力应变图，分析不同深度土层变形特征。2.2 浅部土层本文浅部土层指埋深小于 79.84 m 的土层，包括第一含水层和第一黏土层。其水位变化趋势和变形特征相似，以第一黏土层为例分析（图 3）。该层水位受降水补给和农业开采影响较大（杨会峰等,2021;崔一娇等，2021），水位在一年中随四季波动，27 月降水较少，地下水用水增加，为水位下图2 北京大兴榆垡监测站地层结构及分层标位置示意图Fig.2 Lithology and extensometers at Yufa monitoring station in Daxing,Beijing高校地质学报2 9 卷 4 期594降阶段，8 月至次年 1 月降水增加，为水位回升阶段。第一黏土层在 2008 年至 2012 年间含水层水位持续下降，2012 年至 2016 年受降水增加、开采量减少等因素影响，水位开始抬升；与之对应，土层变形速率在 2012 年出现明显的由快至慢的变化，变形速率由 2012 年之前的 9.54 mm/a 减缓为之后的 4.16 mm/a（图 3a）。从应力应变图（图 3b）中可以看出，2014年之前不管在水位下降阶段还是回升阶段，土层持续压缩，显示出塑性变形特征和变形滞后的效应。变形滞后的原因可能是两方面的，一是土体的蠕变，二是由于弱透水层渗透系数远小于相邻含水层导致的滞后释水。2014 年之后，在水位回升阶段，土层微量回弹，但并不能回弹到压缩之前的位置，存在残余的变形量，在应力应变图中形成回滞环（图 3b），土体既存在弹性变形又存在塑性变形。这说明土体的变形特征随土体压缩变形的发展随时间变化。2.3 中部土层中部土层指埋深在 58.52114.85 m 的土层，包括第二、三含水层和第二黏土层。其水位变化趋势和变形特征相似，以第二黏土层为例分析（图 4）。中部土层水位与浅部土层相似，经历了 2012 年前下降，20122014 缓慢上升，之后呈现下降趋势，但中部土层的年沉降速率无明显变化（图 4a）。第二黏土层年沉降速率基本保持不变，约为 4.81 mm/a。从应力应变图（图 4b）中可以看出，与浅部土层相似，2012 年之前不管在水位下降阶段还图3 第一黏土层变形与水位关系Fig.3 Relationship between the deformation and water level of the first clay layer图4 第二黏土层变形与水位关系Fig.4 The relationship between the deformation and water level of the second clay layer10864220102012201420162018-52468100510 15 20 25 30 354555405060-5051015202530354045505560时间/年水位累计变形累计变形量/mm(a)水位和变形随时间的变化(b)应力应变曲线水位/m水位/m累计变形量/mm141210881012142010201220142016水位/m时间/年累计变形量/mm(a)水位和变形随时间的变化(b)应力应变曲线-5051015202530354045505560657075水位累计变形水位/m累计变形量/mm-5 0 5 10 1520 25 30 354555405060 6575702018595尤斌昊等：北京南部沉降区黏土层变形特征和参数分析4 期是回升阶段，土层持续压缩，显示出塑性变形特征和变形滞后的效应。2012 年之后，在水位回升阶段，土层微量回弹，但不能回弹到压缩之前的位置，存在残余的变形量，在应力应变图中形成回滞环（图 4b），土体既存在弹性变形又存在塑性变形。土体的变形特征随土体压缩变形的发展随时间变化。2.4 深部土层深部土层指埋深大于 114.85 m 的土层，包括第四、五含水层和第三、四黏土层。其水位和变形规律相似，因此只给出第四黏土层水位和变形关系图分析（图 5）。深部土层水位降幅大，第四黏土层水位在 7 年中持续下降，下降速率为 0.65 m/a，土层持续压缩，变形速率 3.90 mm/a（图 5a）。从应力应变图（图 5b）中可以看出，与浅部和中部土层不同，20112017 年间深部土层不管在水位下降还是回升阶段，土层始终持续压缩，没有出现明显的变形回滞环，显示出塑性变形特征和变形滞后的效应。各土层既有弹性变形又存在塑性变形，变形均以塑性变形为主，有明显的变形滞后性。而浅部土层在 2012 年后沉降速率减缓，中部土层保持基本不变的变形速率，浅部和中部土层都显示出弹性变形即可看到土层回弹；深部土层压缩变形速率在 2014 年后呈现加快特征，且始终保持压缩变形，未出现土层回弹阶段。3 黏土层水文地质参数图解法和配线法能很好地利用研究区现有长时间序列水位和沉降观测资料，得出符合野外尺度的水文地质参数，为沉降数学模型的建立和识别提供较合适的初始参数值。下面介绍这两种方法的参数计算过程和结果。3.1 图解法地下水位的下降和上升等效于土层的加载和卸载，土层随之压缩和回弹，形成应力应变曲线的回滞环。图解法就是利用这种回滞环估算相关的水文地质参数（Riley,1969;Burbey,2001）。具体方法是以水位代替有效应力作为为纵坐标，变形代替应变作为横坐标，形成水位变形回滞环，其卸载部分曲线斜率的倒数等于弹性贮水率 Ske，加载部分曲线斜率的倒数等于非弹性贮水率 Skv；过回滞环中最小应力点作斜率为 1/Ske直线，与应变轴有一交点 a，过最小应力点作一条垂直于应变轴的直线，与应变轴有一交点 c，线段 ac的长度即为永久变形量（图 6a）；估算加载过程中的固结度 U=bp/Skv，其中 为应力的变化；选择满足线性压力分布关系的时间因数 Tv，并根据 Kv=SkvbTv/t 算出 Kv。根据北京南部沉降区各黏土层应力应变图，从中截取合适的回滞环（图 6），计算得各黏土层弹性贮水率、非弹性贮水率、垂向渗透系数和变形滞后时间如表 1 所示。因深部第三、四黏土层土层基本不回弹，没有形成回滞环，本次选取的回滞环均来自第一、二黏土层。3.2 配线法配 线 法（Wikramaratna et al.,1985;庄 超 等，图5 第四黏土层变形与水位关系Fig.5 The relationship between the deformation and water level of the fourth clay layer50-5-1020102012201420162018-5-10-505051015202530354540-5051015202530354045时间/年水位累计变形累计变形量/mm(a)水位和变形随时间的变化(b)应力应变曲线水位/m水位/m累计变形量/mm高校地质学报2 9 卷 4 期5962017）假设在一维的含水层系统中：（1）弱透水层是均质的，且渗透系数和贮水率不随时间变化；（2）弱透水层是水平侧向无限延伸的；（3）零时刻（初始时刻）降深在弱透水层中呈线性分布，从零时刻开始下覆含水层降深随时间以速率 线性增大，而上覆含水层水位始终保持不变；（4）含水层系统始终处于饱和状态；（5）弱透水层中的水流为垂向一维流，且服从达西定律。令 td=Kt/(Ssl2)，Fd=2FK/(Ss2l3)，由一维含水层系统的水流控制方程（Wikramaratna et al.,1985）推导可得：F ttntdddd(1 exp(21)=8144n=1(21n)22（1）其中l为弱透水层厚度；K为弱透水层渗透系数；Ss为弱透水层贮水率；为含水层水位降深变化速率。实测的弱透水层累计变形量随时间的变化关系(F-t)曲线和由式（1）所得(Fd-td)关系标准曲线形状相同，只是横纵坐标分别平移了 Kt/(Ssl2)和 2FK/(Ss2l3)。因此，可以运用配线法计算弱透水层的渗透系数 K 和贮水率 Ss。将实测 F-t 曲线匹配 Fd-td标准曲线，在重叠部分取一匹配点可得：K=2Ft lF tdd22（2）Ss=2F tlFtdd（3）以 2011 年作为初始时刻，将各含水层水位降深近似为线性增大，各黏土层实测变形量与标准曲线匹配（图 7），所取匹配点数据及计算渗透系数和贮水率如表 1 所示。3.3 结果图解法结果（表 1）表明第一和第二黏土层非弹性贮水率在 1.4910-52.1010-5之间，弹性贮水率在 3.9310-58.3810-5之间，非弹性贮水率为弹性贮水率的 24 倍；垂向渗透系数在 6.2010-6 4.3910-5之间；土层变形滞后时间为 2.112.8 年。同一土层反演出的变形滞后时间存在一定差异，与图解法中回滞环涉及时间段的应力应变幅度差异有关，通常应力应变幅度大则滞后时间长。配线法估图6 应力应变回滞环（a,b对应第一黏土层中回滞环；c,d对应第二粘土层中回滞环）Fig.6 Stress-deformation loop(a and b correspond to the loops in the first clay layer;c and d correspond to the loops in the second clay layer)(a)(c)(b)(d)597尤斌昊等：北京南部沉降区黏土层变形特征和参数分析4 期算出研究区浅部土层的贮水率在 4.810-4左右，深部土层的贮水率在 1.0810-42.0510-4之间，浅部土层的贮水率是深部土层的 24 倍；浅部土层垂向渗透系数 4.2410-61.0310-5 m/d 比深部土层5.4710-77.7110-7 m/d 大 12 个数量级。对比分析发现，配线法所得贮水率比图解法所得塑性贮水率大 1 个数量级，可能的原因是配线法考虑了黏土层对相邻含水层水位降深变化的滞后响应。由于北京南部沉降区黏土层贮水率、垂向渗透系数、滞后时间等参数未见报道，为检验本参数结果的可靠性、合理性，将其与上海地区弱透水层（岩性为黏土）水文地质参数（叶淑君等，2005）对比（表 1）。可以看出，虽然地区不同，但黏土层的参数值基本在相近的取值范围内。4 结论（1）分析北京南部沉降区榆垡站各土层应力应变图发现，不同深度土层变形特征既有相似之处也存在明显差别。相同之处是各土层既有弹性变形又存在塑性变形，变形均以塑性变形为主，有明显的变形滞后性。区别表现在浅部土层在 2012年后沉降速率减缓，中部土层保持基本不变的变形速率，浅部和中部土层都显示出弹性变形即可看到土层回弹；深部土层压缩变形速率在 2014 年后呈图7 各层累计变形量匹配标准曲线（a,b,c,d分别对应第一到第四黏土层）Fig.7 Observed cumulative compaction data matched with the type curve(a,b,c and d correspond to the first to the fourth clay layers,respectively)表1 黏土层参数值Table 1 Estimated parameter values of clay layers层组图解法配线法Ske Skv Kv/(m/d)(a)SsKv/(m/d)一1.7410-53.9310-56.2010-612.84.810-41.0310-51.4910-54.4710-54.3910-52.1二1.8610-57.3410-51.1110-58.27.4610-44.2410-62.1010-58.3810-59.9510-610.5三-2.0510-47.7110-7四-1.0810-45.4710-7上海地区弱透水层参数经验值1.410-51.410-44.710-49.610-41.110-61.410-410.7241.71.010-67.510-62.910-72.610-4(a)(c)(b)(d)高校地质学报2 9 卷 4 期598现加快特征，且始终保持压缩变形，未出现土层回弹阶段。（2）使用图解法和配线法计算了北京南部沉降区各黏土层弹性贮水率、非弹性贮水率、垂向渗透系数和变形滞后时间等参数。图解法结果表明第一和第二黏土层非弹性贮水率为弹性贮水率的 24倍，土层变形滞后时间为 2.112.8 年。配线法估算出研究区浅部土层的贮水率是深部土层的 24 倍，垂向渗透系数比深部土层大 12 个数量级。两种方法计算所得参数较为接近，计算结果验证了两种方法的可行性，为建立北京南部沉降区区域地面沉降模型提供可参考的水文地质参数值及参数值变化规律。（3）本次研究成果仅依据了榆垡监测站这一个站点的数据，相关结果对北京南部沉降区具有一定参考价值。但要全面揭示南部沉降区土层特征和参数，还需更多分层标监测站数据的支撑。参考文献(Referances):曹群，陈蓓蓓，宫辉力，等.2019.基于SBAS和IPTA技术的京津冀地区地面沉降监测J.南京大学学报(自然科学)，55:381-391.崔亚莉，邵景力，谢振华，等.2003.基于MODFLOW的地面沉降模型研究以北京市区为例J.工程勘察，19-22.崔一娇，杜旋，孙赵爽，等.2021.城市发展模式变化对地下水补给的影响以北京市大兴区为例J.城市地质，16:9-17.董成志，江思珉，夏学敏，等.2017.基于MODFLOW-SUB模块的地层压缩与地面沉降模拟研究J.工程勘察，45:42-48.郭海朋，白晋斌，张有全，等.2017.华北平原典型地段地面沉降演化特征与机理研究J.中国地质，44:1115-1127.郭海朋，李文鹏，王丽亚，等.2021.华北平原地下水位驱动下的地面沉降现状与研究展望J.水文地质工程地质，48:162-171.何庆成，刘文波，李志明.2006a.华北平原地面沉降调查与监测J.高校地质学报，12(2):195-209.何庆成，叶晓滨，李志明，等.2006b.我国地面沉降现状及防治战略设想J.高校地质学报，12(2):161-168.雷坤超，罗勇，陈蓓蓓，等.2016.北京平原区地面沉降分布特征及影响因素J.中国地质，43:2216-2228.李文鹏，王龙凤，郭海朋，等.2021.中国地面沉降防治成效与对策建议J.中国水利，32-35.李文运，崔亚莉，苏晨，等.2012.天津市地下水流-地面沉降耦合模型J.吉林大学学报(地球科学版)，42:805-813.刘妍君，马腾，杜尧，等.2021.黏性土弱透水层压实作用：原理、技术及其水文地质意义J.地学前缘，1-8.刘予，叶超，贾三满.2007.北京市平原地面沉降区含水岩组和可压缩层划分J.城市地质，2:10-15.卢毅，宋泽卓，于军，等.2019.基于BOFDA的砂黏土互层垂向变形物理模型试验研究J.高校地质学报，25:481-486.罗勇，贾三满，赵波，等.2011.北京南部地区地面沉降发育特征及成因分析J.城市地质，6:1-5.罗跃，叶淑君，吴吉春，等.2015.上海市地下水位大幅抬升条件下土层变形特征分析J.高校地质学报，21(2):243-254.罗跃，叶淑君，吴吉春，等.2018.区域地面沉降数值模拟可视化系统开发及应用J.华中科技大学学报(自然科学版)，46:28-33.罗跃，叶淑君，吴吉春，等.2018.地面沉降模型的参数全局敏感性J.浙江大学学报(工学版)，52:2007-2013+2034.孟世豪，崔亚莉，田芳，等.2021.基于MODFLOW-SUB建立变渗透系数的地下水流地面沉降模型研究J.吉林大学学报(地球科学版):1-10.施小清，薛禹群，吴吉春，等.2006.常州地区含水层系统土层压缩变形特征研究J.水文地质工程地质，(33)3:1-6.田芳，郭萌，罗勇，等.2012.北京地面沉降区土体变形特征J.中国地质，39:236-242.田芳，罗勇，周毅，等.2017.北京地面沉降与地下水开采时空演变对比J.南水北调与水利科技，15:163-169.王云龙.2018.中国地面沉降现状图M.北京:地质出版社.熊小锋，罗跃，施小清，等.2017.基于TOUGH2-FLAC(3D)耦合的三维地面沉降数值模拟及控制策略研究J.高校地质学报，23(1):172-180.薛禹群，张云，叶淑君，等.2006.我国地面沉降若干问题研究J.高校地质学报，12(2):153-160.晏霞，刘媛媛，赵振宇.2021.利用时序InSAR技术监测南水进京后北京平原地区的地面沉降J.地球物理学进展，36(6):2351-2361.杨会峰，曹文庚，支传顺，等.2021.近40年以来华北平原地下水位演变及其超采治理建议J.中国地质，48(4):1142-1155.杨勇，李国敏，窦艳兵，等.2010.抽取地下水引起地面沉降的研究现状与进展J.工程勘察，38:32-37+91-93.杨艳，刘贺，罗勇，等.2021.北京东部地区地面沉降发育特征分析J.上海国土资源，42:7-12+18.杨蕴，宋健，朱琳，等.2019.基于KELM地面沉降替代模型的地下水多目标管理模型研究J.南京大学学报(自然科学)，55:349-360.叶超，刘霞，崔文君，等.2018.北京地面沉降M.北京:地质出版社.叶超，田芳，罗勇，等.2019.北京地面沉降控制区划及防控措施J.南京大学学报(自然科学)，55:440-448.叶淑君，薛禹群.2005.应用沉降和水位数据计算上海地区弱透水层的参数J.岩土力学，26(2):256-260.叶淑君，薛禹群，张云，等.2005.上海区域地面沉降模型中土层变形特征研究J.岩土工程学报，27(2):140-147.于军，吴吉春，叶淑君，等.2007.苏锡常地区非线性地面沉降耦合模型研究J.水文地质工程地质，5:11-16.张安京，叶超，李宇.等.2008.北京地下水M.北京:中国大地出版社.张云，薛禹群，吴吉春，等.2006.抽灌水条件下上海砂土层的变形特征和变形参数J.水利学报，37(5):560-566.赵亚丽，王彦兵，王新雨，等.2021.TPCA分析北京平原区地面沉降的时空演化特征J.武汉大学学报(信息科学版)，1-11.周毅，罗郧，郭高轩，等.2016.冲洪积平原地面沉降特征及主控因素以北京平原为例J.地质通报，35:2100-2110.庄超，周志芳，韩江波.2017.一种长期观测资料确定弱透水层参数的方法J.岩土力学，38:1359-1364.庄超，周志芳，李兆峰，等.2017.一种确定超固结弱透水层水力参数的方法J.岩土力学，38:61-66.Burbey T J.2001.Stress-strain analyses for aquifer-system characterization 599尤斌昊等：北京南部沉降区黏土层变形特征和参数分析4 期J.Ground Water,39:128-136.Guo H,Zhang Z,Cheng G,et al.2015.Groundwater-derived land subsidence in the North China Plain J.Environmental Earth Sciences,74:1415-1427.Herrera-Garcia G,Ezquerro P,Tomas R,et al.2021.Mapping the global threat of land subsidence J.Science,371:34-36.Riley F S.1969.Analysis of Borehole Extensometer Data from Central California M/Land Subsidence.Washington:International Association of Hydrological Sciences Publication:423-431.Van Der Kamp G.2001.Methods for determining the in situ hydraulic conductivity of shallow aquitards-an overview J.Hydrogeology Journal,9:5-16