河北省献县地面沉降易发性分区及防治建议

摘要

地面沉降是当今危害较大的一种地质灾害，严重影响区域内建设规划、经济、环境和社会协调发展，对城市发展及人类生命财产安全均有较大危害。献县位于河北省沧州市西部，是地面沉降发育较严重地区，是全国重点蓄滞洪区之一，对献县地面沉降进行易发区划分采取精准防控具有重要意义。本文通过选取深层地下水位、浅层地下水位、第四系厚度、累计沉降量、沉降速率、路网密度6个影响因子，以评价区各乡镇为评价单元，通过灰色关联度法确定地面沉降影响因子权重，利用ArcGIS空间分析功能进行叠加分析，对献县区域进行地面沉降易发性分区，同时提出了防治建议，为地面沉降科学防治提供依据。

关键词

灰色关联度法; 地面沉降; 易发性分区; 影响因子; 权重; 防治建议; 献县

0 引言

地面沉降是由土层压缩而引起地面标高降低的一种常见而缓慢的地质现象^［1］，日益严峻的地面沉降灾害会造成地面高程降低、城市内涝积水、影响南水北调引水工程安全、河道防洪排沥能力降低、影响油田正常开发等一系列灾害^［2-3］，严重制约当地重大工程建设和经济高质量发展。

地面沉降易发区是指地质构造、地层岩性、地下水动态变化和气候等条件容易或者可能引发或加剧地面沉降的区域^［4］。地质灾害易发性评价是以基础地质环境条件（内在控制因素）为出发点，静态考察地质灾害在相对稳定的孕灾环境中发生的可能性大小^［5］。

通过对区域地质环境条件和地质灾害的分布与发育特征的分析，利用模型方法实现地质灾害空间发生的可能性大小的区划，科学评价地质灾害易发性，可为国土规划、地质灾害防治和监测预警提供基础依据^［6］，为合理调整防控措施提供决策依据，为地质灾害防治规划提供科学依据。灰色关联度分析法来源于灰色理论，是根据因素之间发展趋势的相似或相异程度，作为衡量因素间关联程度的一种方法，可以减少人为因素的主观判断。笔者针对献县地面沉降影响因子，运用灰色关联度法建立评价模型，对献县地面沉降发育易发性进行分区，并提出防治建议。

1 区域概况

河北省献县地处冀中平原，位于沧州市西部，东部与沧县相邻，南部与泊头市和衡水市武强县相靠，西部与衡水市饶阳县和肃宁县相接，北部与河间市相连，下辖18个乡镇，总面积1174 km²，全县常住人口56.93万，是京津冀都市圈、环渤海经济圈县之一，经济发展迅速。

献县属冲积扇平原与冲积平原的交接地区，地势平坦，海拔10~15 m，坡降约万分之一。人类工程活动强烈，对地貌形态破坏较大。献县属暖温带半湿润半干旱大陆性季风气候，春季干燥多风，降雨稀少，蒸发强烈，夏季炎热多雨，秋季天高气爽，冬季干燥寒冷。年平均气温11 ℃，年平均降水量560.7 mm。境内主要有滏阳河、滏阳新河、滹沱河、子牙河、子牙新河等5大河流。县境西部滏阳河与滹沱河交汇的三角地带为滞洪区，俗称“四十八村”，是全国重点蓄滞洪区之一。献县全境被第四系覆盖，沉积厚度360~480 m，岩性多为砂、砂砾、粘性土层交错的巨厚松散砂粘土互层。

献县区域按地下水埋藏条件及动力特征，将第四系松散岩类孔隙水分为浅层地下水（潜水或微承压水）与深层地下水（承压水）2种类型。

浅层地下水赋存于全新统（Q₄）地层，底界深度30~40 m，含水层以粉砂为主，少量细砂，厚度7~15 m，条带状分布，单位涌水量<1 m³/h•m，为潜水-微承压水，动态类型为入渗-蒸发型。深层地下水底界埋深480~550 m，含水层岩性为粉砂、细粉砂等，砂层厚度10~60 m，单位涌水量1~5 m³/h•m。其补给来源主要接受侧向径流补给，其次是上覆含水层的越流补给，其动态类型为开采-疏干型。

2 地质灾害发育现状

献县地质灾害发育类型主要为地面沉降，具有影响范围大、持续时间长、发展时间缓慢、成因机制复杂及防治难度大等特点^［7］。根据监测数据，献县地面沉降发育普遍，县域内沉降发育程度不一，沉降速率不同，商林镇累计沉降量已超过1600 mm，沉降量较低区域在1000 mm左右（图1）。近5年献县境内沉降速率最大值发生于2017年，为81 mm/a，近2年整体有减缓趋势，县域内沉降速率发育不均匀，最大沉降速率为50 mm/a以上，最小沉降速率为20 mm/a以下（图2）。

图1  累计地面沉降量

Fig.1  Accumulated land subsidence

图2  沉降速率

Fig.2  Rate of land subsidence

3 地面沉降影响因子

地面沉降影响因子较多且关系复杂，大量研究表明，地下水资源不合理开发利用是诱发地面沉降的主要因素^［8］，第四系厚度是产生地面沉降的内部影响因素，汪洋等^［9］提出第四系地层厚度与地面沉降量大体呈正相关趋势；杜东等^［10］在北京市通州区地面沉降特征与影响因素研究中，提出了路网密度大小与地面沉降严重程度呈正相关。路网密度即区域路网长度与区域面积的比值，用来代表影响地面沉降的动载荷。

充分考虑研究区潜在地面沉降强度和现状地面沉降强度^［11］，结合研究区地质背景、地下水水文地质条件、沉降发展现状及现有社会发展现状，确定研究区6个影响因子，分为潜在地面沉降影响因子与现状地面沉降影响因子^［9］。其中，潜在地面沉降强度影响因子为：浅层地下水水位、深层地下水水位、第四系厚度；现状地面沉降影响因素为：累计沉降量、沉降速率、路网密度。

4 影响因子权重

本文根据灰色关联度原理建立评价模型，确定研究区地面沉降各影响因子权重，减少人为因素主观判断，使权重更客观，定量分析各影响因子间关系^［12-13］ 。

首先确定评价单元，本次评价以献县各乡镇为单位，将评价区分为18个评价单元。将影响因子量化赋值，结合文献资料及献县整体地质情况，按相等组距、相等分裂的方法将影响因子分为4个等级并赋值（见表1）。

表1  地面沉降影响因子及赋值

Table 1  Impact factors and assignments of land subsidence

赋值	影响因子
	累计沉降量/mm	沉降速率/mm	浅层地下水水位埋深/m	深层地下水水位埋深/m	第四系厚度/m	路网密度
1	<1000	<10	<4	<60	<380	<1.5
2	1000~1200	10~30	4~6	60~70	380~420	1.5~2.5
3	1200~1400	30~40	6~8	70~80	420~460	2.5~3.5
4	>1400	>40	>8	>80	>460	>3.5

然后对各影响因子根据2020年沧州市地面沉降监测报告、2020年地下水统测结果、河北省北京市天津市区域地质志和中国分层级道路矢量数据（OSM），将献县各乡镇（评价单元）评分赋值，形成数列，X_i=｛X₀，X₁，X₂，X₃，X₄，X₅｝=｛累计沉降量，沉降速率，浅层地下水水位，深层地下水水位，第四系厚度，路网密度｝，Y_j=｛Y₁，Y₂，Y₃，Y₄，Y₅，Y₆，Y₇，Y₈，Y₉，Y₁₀，Y₁₁，Y₁₂，Y₁₃，Y₁₄，Y₁₅，Y₁₆，Y₁₇，Y₁₈｝=｛张村乡，陌南镇，乐寿镇，商林乡，十五级乡，淮镇，高官乡，临河乡，韩村镇，本斋回族乡，陈庄镇，段村乡，郭庄镇，河城街镇，垒头乡，南河头乡，西城乡，小平王乡｝，由X_i和Y_j形成赋值结果见表2。

表2  地面沉降灾害分值数据

Table 2  Land subsidence disaster score data

评价单元	X0	X1	X2	X3	X4	X5
张村乡	1	4	1	1	3	2
陌南镇	1	4	3	2	4	1
乐寿镇	3	4	1	2	2	2
商林乡	4	4	3	2	2	4
十五级乡	2	3	1	2	1	2
淮镇	1	1	1	2	1	3
高官乡	1	1	1	2	2	2
临河乡	1	4	1	1	4	2
韩村镇	1	2	1	4	2	1
本斋回族乡	2	3	1	2	2	2
陈庄镇	1	4	1	4	2	3
段村乡	2	4	2	2	3	2
郭庄镇	1	2	1	3	1	2
河城街镇	2	3	3	3	1	3
垒头乡	1	2	1	3	2	1
南河头乡	2	4	1	3	2	4
西城乡	3	3	1	2	2	1
小平王乡	1	4	1	1	3	2

将赋值进行无量纲化，即$X_i^{\mathrm{\text{'}}}(k)=X_i(k)/\frac{\mathrm{1}}{n}{\sum }_{k=\mathrm{0}}^n{X}_i(k)$，得到表3。

表3  赋值无量纲化数据

Table 3  Assignment of dimensionless data table

评价单元	X0	X1	X2	X3	X4	X5
张村乡	0.62069	1.28571	0.81818	0.43902	1.38462	0.92308
陌南镇	0.62069	1.28571	1.63636	0.87805	1.84615	0.46154
乐寿镇	1.86207	1.28571	0.81818	0.87805	0.92308	0.92308
商林乡	2.48276	1.28571	1.63636	0.87805	0.92308	1.84615
十五级乡	1.24138	0.96429	0.81818	0.87805	0.46154	0.92308
淮镇	0.62069	0.32143	0.81818	0.87805	0.46154	1.38462
高官乡	0.62069	0.32143	0.81818	0.87805	0.92308	0.92308
临河乡	0.62069	1.28571	1.22727	0.43902	1.84615	0.92308
韩村镇	0.62069	0.64286	0.81818	1.75610	0.92308	0.46154
本斋回族乡	1.24138	0.96429	0.81818	0.87805	0.92308	0.92308
陈庄镇	0.62069	1.28571	0.81818	1.75610	0.92308	1.38462
段村乡	1.24138	1.28571	1.22727	0.87805	1.38462	0.92308
郭庄镇	0.62069	0.64286	0.81818	1.31707	0.46154	0.92308
河城街镇	0.62069	0.96429	1.63636	1.31707	0.46154	1.38462
垒头乡	0.62069	0.64286	0.81818	1.31707	0.92308	0.46154
南河头乡	1.24138	1.28571	0.81818	1.31707	0.92308	1.84615
西城乡	1.86207	0.96429	0.81818	0.87805	0.92308	0.46154
小平王乡	0.62069	1.28571	0.81818	0.43902	1.38462	0.92308

以数列X₀为参比数列，对表3数据计算：

$${\Delta }_{ij}\left(k\right)=│X_i\left(k\right)-X_j(k)│$$,$$(k=\mathrm{1,2},\mathrm{3},\dots \dots ,m)$$

（1）

得到绝对差数列（表4）。

表4  绝对差数列

Table 4  Absolute difference sequence

评价单元	∆01	∆02	∆03	∆04	∆05
张村乡	0.66502	0.19749	0.18167	0.76393	0.30239
陌南镇	0.66502	1.01567	0.25736	1.22546	0.15915
乐寿镇	0.57635	1.04389	0.98402	0.93899	0.93899
商林乡	1.19704	0.84639	1.60471	1.55968	0.63660
十五级乡	0.27709	0.42320	0.36333	0.77984	0.31830
淮镇	0.29926	0.19749	0.25736	0.15915	0.76393
高官乡	0.29926	0.19749	0.25736	0.30239	0.30239
临河乡	0.66502	0.60658	0.18167	1.22546	0.30239
韩村镇	0.02217	0.19749	1.13541	0.30239	0.15915
本斋回族乡	0.27709	0.42320	0.36333	0.31830	0.31830
陈庄镇	0.66502	0.19749	1.13541	0.30239	0.76393
段村乡	0.04433	0.01411	0.36333	0.14324	0.31830
郭庄镇	0.02217	0.19749	0.69638	0.15915	0.30239
河城街镇	0.34360	1.01567	0.69638	0.15915	0.76393
垒头乡	0.02217	0.19749	0.69638	0.30239	0.15915
南河头乡	0.04433	0.42320	0.07569	0.31830	0.60477
西城乡	0.89778	1.04389	0.98402	0.93899	1.40053
小平王乡	0.66502	0.19749	0.18167	0.76393	0.30239

根据表4得到${\Delta }_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{0.01411}$，${\Delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{1.60471}$，计算关联度：

$${\xi }_{\mathrm{i}\mathrm{j}}\left(k\right)=\frac{{\Delta }_{\mathrm{m}\mathrm{i}\mathrm{n}}+\rho {\Delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}}{│X_i\left(k\right)-X_j\left(k\right)│+\rho {\Delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}}$$

（2）

其中，ρ为分辨系数，取值［0，1］，其不影响关联序，一般取值0.5。计算得数列X_j对参考数列X₀的关联系数，见表5。

表5  X_j对X₀关联系数

Table 5  X_j to X₀ correlation coefficient

评价单元	ξ01	ξ02	ξ03	ξ04	ξ05
张村乡	0.55641	0.81659	0.82972	0.52127	0.73905
陌南镇	0.55641	0.44909	0.77045	0.40263	0.84915
乐寿镇	0.59219	0.44223	0.45705	0.46887	0.46887
商林乡	0.40835	0.49520	0.33919	0.34566	0.56740
十五级乡	0.75637	0.66620	0.70041	0.51603	0.72856
淮镇	0.74115	0.81659	0.77045	0.84915	0.52127
高官乡	0.74115	0.81659	0.77045	0.73905	0.73905
临河乡	0.55641	0.57949	0.82972	0.40263	0.73905
韩村镇	0.99022	0.81659	0.42134	0.73905	0.84915
本斋回族乡	0.75637	0.66620	0.70041	0.72856	0.72856
陈庄镇	0.55641	0.81659	0.42134	0.73905	0.52127
段村乡	0.96430	1.00000	0.70041	0.86344	0.72856
郭庄镇	0.99022	0.81659	0.54477	0.84915	0.73905
河城街镇	0.71248	0.44909	0.54477	0.84915	0.52127
垒头乡	0.99022	0.81659	0.54477	0.73905	0.84915
南河头乡	0.96430	0.66620	0.92986	0.72856	0.58023
西城乡	0.48023	0.44223	0.45705	0.46887	0.37063
小平王乡	0.55641	0.81659	0.82972	0.52127	0.73905

根据关联度公式：

$$r_{ij}=\frac{\mathrm{1}}{m}{\displaystyle \sum _{k=\mathrm{1}}^m}{\xi }_{ij}\left(k\right)(\mathrm{3})$$

由表5计算数列X_j对参考数列X₀的关联度，得到表6。

表6  X_j与X₀的关联度

Table 6  The correlation between X_j and X₀

评价单元	r01	r02	r03	r04	r05
关联度	0.71498	0.68826	0.64233	0.63730	0.66552

按式（1）、（2）、（）依次对表3数据以数列X_i为参比数列，求出关联度矩阵R。

$$R=\left[\begin{array}{c}\begin{array}{ccc}\mathrm{1.00000}& \mathrm{0.71498}& \mathrm{0.68826}\\ \mathrm{0.67418}& \mathrm{1.00000}& \mathrm{0.70208}\\ \mathrm{0.62932}& \mathrm{0.68941}& \mathrm{1.00000}\end{array}\begin{array}{ccc}\mathrm{0.64233}& \mathrm{0.63730}& \mathrm{0.66552}\\ \mathrm{0.61702}& \mathrm{0.71988}& \mathrm{0.66269}\\ \mathrm{0.66151}& \mathrm{0.70943}& \mathrm{0.69836}\end{array}\\ \begin{array}{ccc}\mathrm{0.65590}& \mathrm{0.68087}& \mathrm{0.73072}\\ \mathrm{0.62043}& \mathrm{0.74261}& \mathrm{0.74264}\\ \mathrm{0.62526}& \mathrm{0.67079}& \mathrm{0.71441}\end{array}\begin{array}{ccc}\mathrm{1.0000}& \mathrm{0.67476}& \mathrm{0.73509}\\ \mathrm{0.64406}& \mathrm{1.00000}& \mathrm{0.62799}\\ \mathrm{0.68344}& \mathrm{0.60652}& \mathrm{1.00000}\end{array}\end{array}\right]$$

求矩阵R各行平均值${\overline{r}}_{(i)}$，i=0，1，2，……，5。根据式（4）计算各影响因子权重（见表7）。

$$w_i=\frac{{\overline{r}}_{(i)}}{{\sum }_{i=\mathrm{0}}^{\mathrm{5}}{r}_{(i)}}$$

（4）

表7  影响因子权重

Table 7  Influence factor weight

影响因子	X0	X1	X2	X3	X4	X5
权重	0.16554	0.16659	0.16705	0.17045	0.16666	0.16372

根据影响因子权重计算结果，对各因子进行排序，深层地下水位影响权重最大，其次依次为浅层地下水位、第四系厚度、沉降速率、累计沉降量、路网密度。

5 易发区划分

采用综合指数法对研究区进行地面沉降易发分区。

（1）单元栅格划分。将献县县域进行栅格划分，每个栅格大小1 km×1 km，共分为1130个单元栅格。

（2）影响因子指数计算。对各栅格进行影响因子评价赋值，按照综合易发性指数计算公式（5），用加权的方法计算各栅格单元地面沉降易发性指数（G值）。

$$G={\displaystyle \sum _{i=\mathrm{1}}^m}{X}_i·W_i(\mathrm{5})$$

（3）综合指数处理。利用ArcGIS空间分析功能对各影响因素进行叠加分析，得到叠加后研究区的易发性指数和易发指数分布图（图3），经计算得到最高值为3.8911，最低值为1.1737（图4），平均值为2.4455，标准差为0.5195。数值分布较平稳，变异系数较小，说明研究区地质灾害危险性在区域上有一定递变性^［14］，低易发区与高易发区之间有一定过渡区，突变型的区域较少。评价栅格中数值越大，表示地面沉降易发性越大。

图3  地面沉降易发指数栅格图

Fig.3  Grid plot of land subsidence susceptibility index

图4  栅格赋值计算结果分布

Fig.4  Grid assignment calculation result distribution

（4）易发区划分。对易发性指数结合频率和频数，采用自然断点法^［15］，对研究区进行分区，划分为3级，划分区域赋值结果为：高易发区域（1.1737~2.0869），中易发区域（2.0869~2.7230），低易发区域（2.7230~3.8911）。得到易发分区划分和易发分区图（图5）^［16］。

图5  地面沉降易发性分区

Fig.5  Zoning map of land subsidence susceptibility

根据分析结果，献县地面沉降高易发区主要分布于献县西部的陌南镇、临河乡、段村乡、张村乡、商林乡、乐寿镇、小平王乡区域，高易发区面积449.844 km²，占献县总面积的38.40%；中易发区主要分布于献县中部的西城乡、本斋回族乡、十五级乡、南河头乡、陈庄镇区域，中易发区面积421.557 km²，占献县总面积的35.98%；低易发区主要分布于献县东部的韩村镇、垒头乡、郭庄镇、河城街镇、淮镇和高官乡区域，低易发区面积300.156 km²，占献县总面积的25.62%。

6 防控措施建议

深层地下水位是影响地面沉降的首要因素，因此应重点加强对深层地下水开采的管理；尤其对高、中易发区应加强监督，严格控制地下水超采、防止偷采地下水，加大建设供水工程力度；低易发区应加强关注地下水位变化情况，必要时进一步严格地下水资源管理措施^［17］。

同时加强地下水、地面沉降动态监测，对高易发区增加监测点密度，增加监测频次。高易发区应加强城镇市政排水系统建设，防止地面沉降发育带来的城市内涝灾害^［18］。

严格执行地质灾害防治巡排查制度，建立群测群防制度，增强监测预警预报工作。加强工程建设活动监督管理，加强建设项目地质灾害危险性评估，严格审批程序。

7 结论

（1）通过灰色关联度法分析了献县地面沉降各影响因子权重，影响因子中深层地下水位所占权重最大，是影响地面沉降的首要因素。其次是浅层地下水、第四系厚度、沉降速率、累计沉降量、道路密度。

（2）采用综合指数法进行了地面沉降易发分区，将献县划分为高易发区、中易发区、低易发区3个等级。其中，高易发区面积449.844 km²，占献县总面积的38.40%；中易发区面积421.557 km²，占献县总面积的35.98%；低易发区面积300.156 km²，占献县总面积的25.62%。

（3）根据地面沉降易发分区结果，可以指导当地地面沉降防控预案，针对高易发区采取重点防控措施。