贾超, 张少鹏, 丁朋朋, 吴静
山东大学 海洋研究院, 山东 青岛 266237
收稿日期：2020-12-03
基金项目：国家自然科学基金青年基金资助项目(42007234)。
作者简介：贾超(1976-), 男, 江苏徐州人, 山东大学教授。

摘要：针对潍北平原地面沉降问题，建立地面沉降流固耦合三维数值模型，从宏观的角度探讨潍北平原在未来一段时间内区域性地面沉降的发展规律，考虑到地下水开采与地面沉降的相互关系，量化地下水开采对其的负面影响.研究结果表明：潍北平原地区沉降漏斗中心主要分布在寿光市营里镇东部和广陵乡北部；保持2019年地下水开采量，潍北平原未来10年内的沉降仍呈继续发展和增加态势，2029年最大累积地面沉降量将高达565.42 mm，年均增长率为38.97 mm/a，沉降中心未发生转移；随着地下水开采量的减少，地面沉降发育程度明显得到缓解，随之持续减小，地面沉降的减缓趋势则会逐渐趋于平缓；当减采量控制在10%与30%之间时，地面沉降灾害的减缓效果较佳.
关键词：潍北平原地面沉降地下水开采数值模型演化趋势
Temporal-Spatial Response of Land Subsidence to Groundwater and Its Non-engineering Subsidence Control
JIA Chao, ZHANG Shao-peng, DING Peng-peng, WU Jing
Institute of Marine Science and Technology, Shandong University, Qingdao 266237, China
Corresponding author: JIA Chao, E-mail: jiachao@sdu.edu.cn.

Abstract: Aiming at the problem of land subsidence in Weibei Plain, a three-dimensional numerical model of fluid-solid coupling for Weibei Plain land subsidence was established. The development law of Weibei Plain land subsidence in the future was discussed from a macro perspective, and the negative impact of groundwater exploitation on the ground subsidence was quantified by taking into account the mutual relationship between groundwater exploitation and land subsidence. The results show that the subsidence funnel center is mainly distributed in the east of Yingli Town and the north of Guangling Township. If the amount of groundwater exploitation in 2019 is maintained, the subsidence of Weibei Plain will continue to develop and increase in the next 10 years, and the maximum cumulative land subsidence will be as high as 565.42 mm in 2029, with an average annual growth rate of 38.97 mm/a, and the subsidence center will remain unchanged. With the decrease of groundwater exploitation, the development degree of land subsidence can be alleviated obviously, and will continue to decrease, causing the mitigation trend of land subsidence. When the reduction is controlled between 10% and 30%, the mitigation effect of land subsidence disaster could be better.
Key words: Weibei Plainland subsidencegroundwater exploitationnumerical modelevolution trend
随着工业的快速发展及人类社会需求的增长，过度开采地下水诱发的地面沉降已经成为了一个使地表高程降低的全球性环境地质灾害^[1].实践经验表明，地下水超采是造成地面沉降的主要原因^[2-3].
墨西哥城是典型的地面沉降灾害地区，自1891年首次受到关注以来^[4]，全球范围内许多国家，特别是人口密集的地区，均有不同程度的地面沉降灾害记录^[5-7].
潍北平原位于山东半岛蓝色经济区，拥有中国较大规模的蔬菜种植和卤水化工基地，地下水长期处于超量开采状态，地面沉降问题日趋突出^[8].缓解地下水开采引起的地面沉降灾害措施主要分为减少地下水开采量、人工回灌和禁采三种方法^[9-10].最有效且经济的非工程遏制措施就是减少地下水尤其是承压水的开采量.
关于地面沉降的研究目前相对较侧重于对地面沉降模拟方式的改进，或者是探讨土体本构关系及参数变化^[11-12]对地面沉降的影响特征和规律.事实上，从宏观的角度研究大区域尺度下的地面沉降对地下水开采的响应程度，对于区域性地面沉降防范与治理具有十分重要的意义.
为更好掌握潍北平原地面沉降发展趋势，本文以潍北平原为研究背景，根据研究区详细的地质资料，结合高精度的地下水位和地面沉降监测数据，系统分析潍北平原地面沉降分布特征；以达西定律和Biot固结理论为基础，建立潍北平原地面沉降多孔介质三维流固耦合数值模型，模拟预测未来10年内的地下水与区域性地面沉降的时空动态演化过程；在保证地下水开采价值的前提下，探讨地下水减采在何种程度上可以有效降低沉降的负面影响，提出合适的非工程控沉措施.
1 理论与方法地面沉降就是大量抽取地下水改变渗流场从而引起土体有效应力增加，应力场的改变引起位移场发生相应变化.同时需要注意的是，应力和位移的改变将引起土体孔隙结构改变，从而反映出土体孔隙率和渗透系数的变化，导致渗流场再次发生变化，如此反复，地面沉降计算就是地下水渗流模型和土体变形模型相结合的结果.
1.1 地下水三维渗流数学模型考虑达西定律和质量守恒的连续性方程，取多孔介质各向异性主渗透方向与坐标轴方向一致，建立多孔介质三维非稳定地下水渗流模型^[13]，即

(1)

式中：k_x，k_y，k_z分别为各主方向渗透系数；h为点(x，y，z)在t时刻的水头值；W为源汇项；μ为储水率；t为时间.
1.2 地面沉降流固耦合数学模型Biot^[14]在1941年依据严格的固结机理推导出能够准确反映孔隙水压力消散与土体骨架变形相关的“真三维固结理论”.假设土骨架变形符合线弹性理论，变形微小，土颗粒和孔隙水均不可压缩，渗流条件符合达西定律，根据有效应力原理，总应力为有效应力σ′与孔隙水压力u之和，孔隙水不受剪应力，用位移和孔隙压力表示的土体平衡微分方程为

(2)

式中：G为剪切模量；ν为泊松比；γ为土的重度；为拉普拉斯算子；W_x，W_y，W_z分别表示土体在x，y，z方向上的位移量.
根据假设，土体为完全饱和，土粒和水都属于不可压缩，渗流符合达西定律，对于土体微分单元体d_x d_y d_z，某一时间段内流出的水量等于该时间段内单元体的体积变化，则用位移和孔隙水压力表示的渗流连续性方程为

(3)

式中，γ_w为水的容重.
在地面沉降过程中，土体的变形会对孔隙率和渗透系数产生一定的影响，根据Kozeny-Carman方程，其动态变化过程可描述为

(4)

式中：n₀为土体初始孔隙率；n为土体变形后孔隙率；k₀为土体初始渗透系数；k为土体变形后渗透系数；ε_v为土体的体应变量，为线应变之和.
将式(2) 和式(3) 联立建立Biot三维固结方程组.结合式(4)、初始条件和边界条件^[15]就可以建立地下水开采与地面沉降耦合数值模型进行求解.
初始条件主要包括位移初始条件式(5)、孔隙水压力初始条件式(6) 和地应力初始条件式(7)，即

(5)

(6)

(7)

式中：u₀(x, y, z) 为已知初始孔隙水压力；σ_x，σ_y为土体的初始水平应力；σ_z为土体的初始垂向应力；z为计算点深度；K₀为静止侧压力系数.
地下水边界条件一般包括三类：Dirichlet边界、Neumann边界和混合边界，其表达式分别见式(8) ~式(10)：

(8)

式中，H₁(x, y, z)为水头边界Γ₁上的已知水头函数.

(9)

式中，q(x, y, z, t) 为水头边界Γ₂上的已知流量.

(10)

式中：α，β为边界Γ₃上的已知数；H为边界Γ₃上的自由水头.
对于开采井边界，使用单位长度的质量流量N₀表示：

(11)

式中：N₀为线源质量流量，负值表示流出水量；κ为土体的渗透率；μ为水的动力黏度；ρ_w为水的密度；U₀为开采井水量流出速度；▽_u为孔隙水压力u的哈密顿算子.
2 算例分析2.1 研究背景潍北平原涉及地区主要为山东省潍坊市的寿光市、寒亭区、昌邑市，总面积约为4 096 km²(图 1).地形平坦微向北倾斜，坡降为1‰~3‰，南部有残丘零星分布，西部与华北平原接壤，南部为低山丘陵区并逐渐向鲁中山地过渡，北部为河流冲洪积及海积平原区，揭露的土体岩性主要为第四系冲洪积形成的中粗砂、砾石以及冲海积形成的粉细砂、亚砂土、亚黏土.
图 1(Fig. 1)

图 1 潍北平原地区及监测点地理位置图Fig.1 Geographical location of the Weibei Plain and the monitoring points (a)—潍北平原地理位置; (b)—沉降、地下水监测点标识; (c)—图 1a中研究区交通、监测点、边界.

2010年至2018年平均降水量为577.67 mm，平均蒸发量为1 834.0 mm.潍北平原研究区北部为莱州湾南岸海岸线，西部分布有小清河及塌河，东部有胶莱河、弥河、白浪河等贯穿内部.
研究区南部山前冲洪积平原主要开采浅层含水层内的淡水，用于农业灌溉和城镇生活用水；北部滨海海积平原主要开采中深层地下水用于卤水工业，居民生活主要饮用地表水.
2.2 地面沉降分布特征基于InSAR技术获取的潍北平原地面沉降监测数据，按沉降速率将研究区内沉降区域划分为不同发育程度等级，潍北平原2017—2019年地面沉降速率分区情况如图 2所示.潍北平原地面沉降强发育区主要分布在寿光市和昌邑市，沉降速率大于50 mm/a的区域面积约为302.04 km²，沉降速率30~50 mm/a的区域面积约为1 308.64 km².
图 2(Fig. 2)

图 2 潍北平原2017—2019年地面沉降速率分区图Fig.2 Zoning map of land subsidence rate in Weibei Plain from 2017 to 2019

为准确监测分析各区域的地面沉降发展演化规律，山东省某单位于2012年在山东省重点沉降区域布设了地面沉降监测系统，其在潍北平原共有9个监测点(图 1c)，其中JC1，JC2和JC3监测点分别在广陵乡以北约6.3 km处、营里镇以北约3.2 km处和羊口镇以南约5.5 km处，均位于潍北平原北部沉降严重区域，并分别于2017年、2018年和2019年在该地区进行了3次二等水准测量，2017年至2019年间水准测量点的年平均沉降速率最大为92 mm·a^-1.
据调查研究，潍北平原中深层地下水降落漏斗与地面沉降分布特征之间呈现出较好的相关性^[16]，地面沉降漏斗区域与中深层地下水降落漏斗区域基本保持一致，由此推断潍北平原地区产生地面沉降的主要原因是由于中深层地下水的过量开采.
2.3 三维流固耦合数值模型以达西定律和Biot固结理论为基础，在计算得出模型水力边界条件和抽水井水力边界条件的前提下，通过求解耦合多孔弹性介质骨架的地面沉降三维流固耦合数值模型，可以求得模型的变形和应力，从而计算得到地面沉降变形量.
2.3.1 水文地质概念模型对于浅层含水层组(Ⅰ，Ⅱ，Ⅲ层)，将北部莱州湾南岸海岸线设为定水头边界，西部小清河与塌河设为定流量边界，东部胶莱河设为河流边界，南部由于接收山前侧向径流补给，设为定水头边界；对于中深层含水层组(Ⅳ，Ⅴ，Ⅵ层)，其与浅层含水层组无明显水力联系，南部设为定流量边界，其余边界条件假设与浅层含水层组一致.垂向上表面要经过大气降水补给、河流侧渗补给及蒸发，是水量不断交换的边界；下部为黏土层，设为隔水边界.
根据现有的水文地质资料，在满足计算精度的前提下，将模型分成639 189个的立方体网格，如图 3所示.
图 3(Fig. 3)

图 3 水文地质概念模型Fig.3 The hydrogeological conceptual model

2.3.2 流固耦合数值模型根据已掌握的地层资料和当地抽水试验数据，以及本次研究区地下水及矿产开采相关资料，建立三维几何数值模型(图 4).由于开采井径相对于整个研究区面积十分微小，为了在不影响精度的前提下提高计算效率，将对应的开采流量换算成线源质量流量的形式作为边界条件，将372口不同地理位置的开采井及不同开采条件赋值到模型中，表 1列举了研究区部分开采井地理位置、年均质量流量，其中土体参数如表 2所示.
图 4(Fig. 4)

图 4 三维流固耦合数值模型Fig.4 The three-dimensional fluid-solid coupled numerical model

表 1(Table 1)

表 1 研究区部分开采井位置及开采量Table 1 Location and the corresponding flow of part of the wells 编号 坐标 年均质量流量N0 kg·m-1·s-1 k1 118.879 9, 37.003 6 0.030 7 k2 118.824 3, 37.055 9 0.006 6 k3 118.952 7, 37.057 6 0.014 9 k4 118.749 1, 36.974 3 0.126 1 k5 118.918 7, 36.924 3 0.005 8 k6 118.818 1, 37.106 9 0.088 8 k7 118.926 3, 36.858 4 0.043 1 k8 118.919 7, 36.919 1 0.065 5 k9 118.959 6, 37.025 0 0.095 4 k10 119.112 0, 37.118 1 0.008 3 k11 119.394 6, 37.027 7 0.054 8

表 1 研究区部分开采井位置及开采量 Table 1 Location and the corresponding flow of part of the wells

表 2(Table 2)

表 2 潍北平原地区地层物理力学参数Table 2 Physical and mechanical parameters of the strata in the Weibei Plain area 地层 底板埋深d/m 弹性模量E/MPa 泊松比 密度ρ/(kg·m-3) 孔隙率 压缩系数α1-2/MPa-1 渗透系数k/(m·d-1) Ⅰ 6~10 49~52 0.37~0.39 1 985~1 998 0.412~0.439 — 0.92~1.23 Ⅱ 27~39 51~56 0.36~0.38 2 002~2 026 0.392~0.408 0.1~0.4 0.007~0.009 Ⅲ 43~61 77~83 0.23~0.25 2 025~2 039 0.376~0.388 — 3.22~3.47 Ⅳ 88~93 59~62 0.29~0.32 2 041~2 042 0.377~0.385 0.18~0.38 0.005~0.008 Ⅴ 116~120 120~125 0.19~0.22 2 047~2 054 0.372~0.382 — 1.83~2.31 Ⅵ 149~152 56~62 0.28~0.31 2 069~2 071 0.355~0.365 0.2~0.35 0.002~0.003

表 2 潍北平原地区地层物理力学参数 Table 2 Physical and mechanical parameters of the strata in the Weibei Plain area

北部的莱州湾南海岸线、东部的胶莱河设为定水头边界；南部界限和西部的小清河及塌河设为流量边界，其中流量边界水力变化条件通过地下水流场动态演化分析计算得到，同时限制水平方向移动，竖直方向自由滑动；模型上表面为自由边界，底部为隔水边界，水平和竖直方向固定位移；开采井采用流量边界条件，以实际开采量转化为单位质量流量，限制水平方向位移，竖直方向自由移动.
2.4 模型校核与验证2.4.1 地下水数值模型校核将2016年6月至2018年9月作为地下水模型识别期，2018年10月至2019年10月作为模型检验期.各水位监测点处的浅层和中深层地下水水位实测值与计算值拟合验证结果见图 5，其R²分别为0.995，0.997，计算值拟合效果良好.
图 5(Fig. 5)

图 5 地下水位散点图Fig.5 Scatter plot of the groundwater level (a)—浅层地下水监测点; (b)—中深层地下水监测点.

选取寿光市广陵乡北部沉降漏斗处的Q27号浅层地下水和S26号中深层地下水监测井(图 1b，图 1c) 处实测数据与模拟计算值进行水位动态校核，如图 6所示.其最大差值分别为0.387 m和0.357 m, 拟合误差均小于拟合计算期间内水位变化值的10%.
图 6(Fig. 6)

图 6 地下水位动态校核图Fig.6 Dynamic check of the groundwater level (a)—Q27号浅层地下水监测井; (b)—S26号中深层地下水监测井.

检验期内边界流入总量为18 681.541 5万m³/d，流出总量为18 681.723 8万m³/d，差值为-1 821.98 m³/d，误差为-0.000 975%，地下水量呈均衡状态.该地下水模型具有一定的准确性和可靠性.
2.4.2 地面沉降趋势校核通过将模拟沉降结果与依据GIS空间分析技术处理得到的2017—2019年地面沉降InSAR监测数据进行对比分析(图 7)，确定三维流固耦合数值模拟沉降在宏观上时空分布趋势的准确性.
图 7(Fig. 7)

图 7 潍北平原2017—2019年累积地面沉降量分布云图Fig.7 Cloud picture of the cumulative land subsidence in the Weibei Plain from 2017 to 2019 (a)—InSAR监测数据; (b)—数值模拟结果.

根据图 7分析研究区地面沉降模拟结果分布特征，沉降区域主要发生在寿光市广陵乡北部、营里镇东部、羊口镇、寒亭区大家洼镇、央子街道，昌邑市下营镇.其中寿光市广陵乡北部地区处的地面沉降漏斗面积最大，影响最严重.从数值模拟结果和InSAR监测数据对比可以看出，在该时间段的地面沉降模拟效果较为理想，能够反映此段时间内该地区的实际沉降趋势.
2.4.3 地面沉降水准监测点校准潍北平原研究区内目前设有9个有效二等水准监测点(图 1b，图 1c)，其基本位于沉降较严重的区域.将各水准监测点2017—2019年实测年平均沉降速率与数值模拟结果进行统计对比分析，定量校准潍北平原地面沉降模型的准确性和有效性(图 8).
图 8(Fig. 8)

图 8 水准监测点监测值与模拟值Fig.8 Measured values and simulated values of the monitoring points

对比各监测点实际监测值和模拟值，可以发现模拟值和监测值较符合，最大的模拟偏差为7.8%，最小的模拟偏差为-0.94%，沉降模拟误差均在10%以内.其中位于潍北平原北部沉降严重区域的JC1，JC2和JC3监测点，模拟年均沉降速率分别达到了86.3，72.3和49.5 mm·a^-1，其沉降速率较明显，表明该地正以较快的速度形成区域性地面沉降大漏斗.
综合各监测点沉降量及地面沉降演化趋势与模拟结果比较分析，认为所建立潍北平原地面沉降三维流固耦合数值模型满足精度要求，整体拟合效果较好，能够反映该地区的地面沉降发展演化趋势和规律.
2.5 地面沉降趋势预测为有效防控潍北平原地面沉降趋势，提出合理的地下水开采方案，运用已建立的潍北平原地面沉降三维流固耦合数值模型，设计现状开采条件、减小地下水开采量等开采方案，预测分析不同地下水开采设计方案对地面沉降特征及规律的影响.
2.5.1 现状开采地下水假设现状地下水开采的条件不变，选取2019年地面沉降量作为初始状态，将2019年的地下水开采量作为现状开采量，保持模型位移边界条件不变，水力边界条件根据地下水流场模拟预测结果施加，预测研究区未来10年内地面沉降发展规律.其中2019年，2025年，2029年的潍北平原地面沉降预测如图 9所示.
图 9(Fig. 9)

图 9 现状开采条件下地面沉降量模拟预测Fig.9 Prediction of land subsidence under current groundwater exploitation condition (a)—2019年; (b)—2025年; (c)—2029年.

由图 9分析可知，如果保持目前开采状态，未来10年内，地面沉降受灾区位置基本未发生明显变化，沉降仍呈继续发展和增加态势，寿光市广陵乡北部(P1) 和营里镇东部(P2) 两处较大的地面沉降漏斗都会相互影响，逐渐相连并向四周扩大.2029年时，研究区最大累积地面沉降量将会达到565.42 mm，相比2019年增加了389.79 mm，年均增长率为38.97 mm·a^-1.根据开采井分布资料，地面沉降基本发生在北部开采井集中分布和开采强度大的区域，而研究区北部开采井主要开采中深层地下水，说明潍北平原地面沉降的发生是开采中深层地下水造成的，大量开采地下水势必会导致并加剧地面沉降的发生.
如图 10所示，潍北平原地面沉降主要发生在寿光市古城街道、营里镇东部、广陵乡北部、侯镇西部、羊口镇南部、寒亭区央子街道东南部、固堤镇北部和昌邑市下营镇.此外，在寿光市城区、稻田镇北部、台头镇也存在不同程度的沉降，其中寿光市广陵乡地面沉降漏斗为影响面积最大、沉降最严重的单个漏斗.到2029年时，寿光市营里镇和广陵乡处的地面沉降漏斗已经构成一个大型区域性漏斗，其影响面积达到953.767 km².
图 10(Fig. 10)

图 10 现状开采下2029年地面沉降漏斗位置图Fig.10 Location map of land subsidence funnel in 2029 under current groundwater exploitation condition

2.5.2 减采地下水以2019年各开采井地下水开采量为基础，设计开采量分别减小10%，30%和50%三种地下水开采方案，通过潍北平原地面沉降流-固耦合模型进行模拟计算得到研究区2029年时的地面沉降分布图，以期分析不同地下水控采量对地面沉降的响应，提出适宜的非工程控制沉降措施.
图 11为潍北平原在2019年地下水开采量的基础上分别减小10%，30%和50%情况下2029年时的累积地面沉降量分布演化云图.结合图 9分析，地下水现状开采、减采10%、减采30%及减采50%等四种工况下的研究区最大累积地面沉降量分别为565.42，518.28，423.44和328.61 mm，各减采工况下其最大累积地面沉降量降幅分别为8.34%，25.11%和41.88%，可见研究区累积地面沉降量较未采取地下水减采措施时明显降低，其地面沉降发育程度得到较大的缓解，与2019年时研究区地面沉降发展形式相比，虽地面沉降影响区域有所扩大，依旧围绕P1与P2地区发展，但程度有所减缓.
图 11(Fig. 11)

图 11 减采作用下2029年地面沉降量模拟预测图Fig.11 Prediction of land subsidence in 2029 under reduction of groundwater exploitation (a)—减采10%; (b)—减采30%; (c)—减采50%.

特别地，在现状开采条件下2025年时研究区最大累积地面沉降量为428.91 mm，超过减采30%工况下2029年沉降值5.47 mm，说明减采地下水这一措施能够有效地减缓地面沉降发生的程度，若不及时采取合理的方案，则地面沉降灾害将会难以抑制.
为做好针对性控沉措施，对不同地下水开采工况下潍北平原地区2029年累积地面沉降量超过300 mm及以上的区域(地面沉降中等发育区) 面积(图 12) 进行提取并分析研究.
图 12(Fig. 12)

图 12 2029年沉降区域面积统计图Fig.12 Statistics of subsidence area in 2029

地下水现状开采、减采10%、减采30%和减采50%四种工况下地面沉降中等发育区面积随开采量的减小而显著减小，其值分别为298.65，247.69，141.31和48.31 km²，面积占比分别为7.38%，6.12%，3.49%和1.19%.减小地下水开采能够有效抑制地面沉降严重区扩展，但过小的减采量对地面沉降的防治效果又会减弱，当减采量超过30%左右时，相比持续减采时面积占比仅减小了2.3%.
统计图 10中A-A′剖面上2029年时现状开采、减采10%、减采30%和减采50%四种工况下的累积地面沉降量值，见图 13，该断面包含了研究区最严重的寿光市广陵乡地面沉降漏斗.从图 13中可以明显看出，保持现状开采地下水量时，潍北平原北部地区地面沉降量最大为555.28 mm，南部沉降量最大为54.98 mm，北部地面沉降发展情况相对于南部较严重，地面沉降漏斗中心主要集中在寿光市广陵乡北部.
图 13(Fig. 13)

图 13 A-A′剖面2029年累积地面沉降量Fig.13 Cumulative land subsidence of section A-A′ in 2029

随着地下水开采量的减少，沉降得到明显缓解，但地面沉降漏斗中心位置依旧不变.对于地面沉降发生较轻的南部地区，例如稻田镇地区，减小开采量这一措施并没有较明显的成效，特别是减采10%地下水开采量措施对累积沉降量小于100 mm的地区几乎没有沉降变化的影响.
通过分析对比各工况下对地面沉降的缓解规律，可以发现地下水控采量并不是越大或者越小就越好.随着开采量的减小，累积地面沉降量均会显著减小，但随着开采量的继续减小，累积地面沉降量的减小趋势会逐渐趋于平缓.提取沉降漏斗中心处各工况下的沉降数据分析，现状开采至减采10%地下水、减采10%至30%地下水和减采30%至50%地下水工况下地面沉降量差值分别为49.48，94.80和72.49 mm，其沉降减缓速率并不随开采的增加而增加，减采30%至50%时沉降变化值较减采10%至30%时小22.31 mm，说明当减采量控制在10%~30%之间时，其减缓地面沉降灾害效果较佳.
潍北平原处于山东半岛蓝色经济区，对地下水资源有较大的需求量，地面沉降防治的同时需要考虑当地经济的可持续发展和人类生活的物质需求.为了能更好地减缓甚至防范地面沉降的发生，本文认为在2019年地下水开采量的基础上减采10%~30%是合理有效的地下水开采方案.同时，地下水的减采应是“开源节流”，对于部分地区应该具有针对性.当地的决策者需要制定严格的环境保护条例和地下水开采方案，公民们也应该提高保护水资源的意识，不能浪费地下水资源，共同营造人与自然和谐共生的氛围.
3 结论1) 保持2019年地下水开采量不变，潍北平原未来10年内的地面沉降灾害仍呈继续发展和增加的态势，区域性地面沉降漏斗分别以寿光市营里镇东部和寿光市广陵乡北部为中心逐渐变大并相连.
2) 当地下水减采量增大时，研究区地面沉降发育较未采取措施时明显得到遏制，但沉降中心并不随地下水开采量的减少而改变地理位置.
3) 地面沉降并不随地下水开采量的减小而线性减缓.随着开采量的减小，沉降漏斗处累积地面沉降量均显著减小，但随着开采量的继续减小，累积地面沉降量的减小趋势逐渐趋于平缓，对沉降发育程度中等及以上面积占比进行分析也能得到相同的规律.
4) 地下水减采能够有效地控制当地地面沉降的发展趋势，当减采量控制在10%~30%时，其减缓地面沉降灾害效果较佳.同时需要实施针对性的地下水阶梯式减采方案来实现地面沉降预防治理.
参考文献

[1]	Galloway D L, Erkens G, Kuniansky E L, et al. Preface: land subsidence processes[J]. Hydrogeology Journal, 2016, 24(3): 547-550. DOI:10.1007/s10040-016-1386-y
[2]	Xue Y Q, Zhang Y, Ye S J, et al. Land subsidence in China[J]. Environmental Geology, 2005, 48(6): 713-720. DOI:10.1007/s00254-005-0010-6
[3]	周志芳, 郑虎, 庄超. 论地下水资源的永久性消耗量[J]. 水利学报, 2014, 45(12): 1458-1463. (Zhou Zhi-fang, Zheng Hu, Zhuang Chao. Study on the unrecoverable depletion of groundwater resource[J]. Journal of Hydraulic Engineering, 2014, 45(12): 1458-1463.)
[4]	Ortega G A, Rudolph D L, Cherry J A. Analysis of long-term land subsidence near Mexico City: field investigations and predictive modeling[J]. Water Resources Research, 1999, 35(11): 3327-3341. DOI:10.1029/1999WR900148
[5]	Rateb A, Zaki A. Inferencing the land subsidence in the Nile Delta using Sentinel-1 satellites and GPS between 2015 and 2019[J/OL]. Science of the Total Environment, 2020, 729(2020-08-10)[2020-10-01]. https://www.sciencedirect.com/science/article/abs/pii/S0048969720323858.
[6]	Chen D, Lu Y, Jia D. Land deformation associated with exploitation of groundwater in Changzhou City measured by COSMO-SkyMed and Sentinel-1A SAR data[J]. Open Geosciences, 2018, 10(1): 678-687. DOI:10.1515/geo-2018-0054
[7]	Jia C, Zhang Y W, Han J M, et al. Susceptibility area regionalization of land subsidence based on extenics theory[J]. Cluster Computing, 2017, 20(1): 1-14.
[8]	Xue Y Q, Wu J C, Xie C H, et al. Research of seawater and saltwater intrusion of Laizhou Bay[J]. Chinese Science Bulletin, 1997, 11(22): 2360-2368.
[9]	Galloway D L, Burbey T J. Review: regional land subsidence accompanying groundwater extraction[J]. Hydrogeology Journal, 2011, 19(8): 1459-1486. DOI:10.1007/s10040-011-0775-5
[10]	Ali M Z, Chu H J, Burbey T J. Mapping and predicting subsidence from spatio-temporal regression models of groundwater-drawdown and subsidence observations[J]. Hydrogeology Journal, 2020, 28(8): 2865-2876. DOI:10.1007/s10040-020-02211-0
[11]	Li Z, Luo Z J, Wang Q, et al. A three-dimensional fluid-solid model, coupling high-rise building load and groundwater abstraction, for prediction of regional land subsidence[J]. Hydrogeology Journal, 2019, 27(4): 1515-1526. DOI:10.1007/s10040-018-01920-x
[12]	Tsai T L, Jang W S. Deformation effects of porosity variation on soil consolidation caused by groundwater table decline[J]. Environmental Earth Sciences, 2014, 72(3): 829-838. DOI:10.1007%2Fs12665-013-3006-7
[13]	Bathe K, Khoshgoftaar M R. Finite element free surface seepage analysis without mesh iteration[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1979, 3(1): 13-22.
[14]	Biot M A. General theory of three-dimensional consolidation[J]. Journal of Applied Physics, 1941, 12(2): 155-164.
[15]	唐金颖. 地下水开采引起地面沉降的数值模拟分析[D]. 沈阳: 东北大学, 2008. (Tang Jin-ying. Numerical simulation and analysis of land subsidence by groundwater exploitation[D]. Shenyang: Northeastern University, 2008. )
[16]	丁朋朋, 贾超, 韩瑶, 等. 时序InSAR技术的潍北平原地面沉降分析[J/OL]. 测绘科学, 2021, 46(2020-09-10)[2020-10-05]. https://kns.cnki.net/kcms/detail/11.4415.P.20200908.1921.002.html. (Ding Peng-peng, Jia Chao, Han Yao, et al. Analysis of land subsidence in Weibei Plain based on time series InSAR method[J/OL]. Science of Surveying and Mapping, 2021, 46(2020-09-10)[2020-10-05]. https://kns.cnki.net/kcms/detail/11.4415.P.20200908.1921.002.html. )