近40年三江平原典型区土地利用变化及其对地下水的影响

束龙仓，王 哲，袁亚杰，张凤海，刘 武，鲁程鹏

（1.河海大学水文水资源学院，江苏 南京 210098；2.河海大学水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；3.黑龙江省水文水资源中心，黑龙江 哈尔滨 150001）

三江平原作为我国重要的商品粮基地和粮食战略储备基地，在我国农业发展中占有举足轻重的地位，其灌溉水源主要采用地下水，而近些年来随着水稻灌区生产规模的不断扩大，过度依赖地下水灌溉使得建三江垦区形成了大型地下水降落漏斗雏形，地下水环境的快速改变严重制约了三江平原灌区农业的可持续发展。同时，三江平原拥有多处湿地自然保护区，地下水资源也是维持保护区的重要因素之一［1］。为了解决该地区发展所面临的地下水水位持续下降的难题，保障当地地下水资源的可持续利用，支撑该地区经济社会的可持续发展，亟待对三江平原地下水动态变化规律进行相关研究，而计算水均衡是揭示地下水动态变化规律的关键。

近年来，国内外相关学者围绕水均衡要素进行了较为深入的研究，探究了如河道硬化［2］、气候变化［3-5］及季节性冻土［6］等对水均衡的影响；同时还在研究中采用White法［7］估算地下水蒸散发量；通过建立地表水—地下水分布式耦合模型［8］、采用秩次集对法［9］模拟和预测地表水与地下水动态变化；利用分布式模拟模型［10］探究地下水在“自然—社会”二元水循环过程中的补给规律等，这些研究为流域水均衡提供了关键的研究思路和方法。虽然诸多学者［11-14］采用不同方法对水均衡要素进行了分析与研究，但往往忽略了研究区下垫面条件的改变对水均衡的影响，土地变化/覆被变化（Land use/cover change，LUCC）是人类活动最直接的表现形式，作为人地关系研究的重要组成与核心内容，也是当前国内外研究热点问题之一［15］，却很少有学者研究其对地下水位埋深和流域水均衡状态的影响。

本文选取三江平原典型区1980、1990、2000、2010及2018年五期的遥感影像数据进行土地利用动态分析与时空迁移分析，揭示近40年流域土地利用的变化特征，探究土地利用类型变化对地下水位埋深的影响以及城市化进程与耕地面积变化对流域主要水均衡要素的影响，为保障该地区的粮食用水安全提供科学依据，为三江平原地下水的可持续开发利用提供理论支撑。

1 研究区与数据

1.1 研究区概况 三江平原典型区位于黑龙江省东北部，是由松花江、黑龙江、乌苏里江以及挠力河所围成的建三江区域，包括富锦市、同江市、抚远市以及饶河县挠力河以北部分，地理坐标范围为44°48′—48°23′N，131°26′—134°47′E，总面积约2.21万km2（图1）。典型区属寒温带湿润-半湿润大陆性季风气候，多年平均降水量为383.5～886.1 mm，降水年内分布不均，主要分布在6—9月，约占全年降水量的70%，全区多年平均气温为2.8 ℃。

典型区内水系发育，河流纵横，均属黑龙江水系，其主要支流有松花江、乌苏里江，区内湖泡星罗棋布，沼泽湿地发育，有大小河流100余条，均以黑龙江、松花江、乌苏里江三个水系为主干。

典型区含水层结构类型主要为前第四纪基岩单层结构含水层系统和东部黏土质低平原双层结构孔隙含水层系统，其中第四系以河流冲积物为主，地表普遍被粉质黏土、黏性土及砂砾质土覆盖。松散岩类孔隙水是区内分布范围最广泛的地下水类型，基岩裂隙水分布范围较小，富水性差。区内富锦集贤以东的广大地区地下水水量丰富，含水层岩性主要为粉质黏土、粉质砂土及细砂，含水层厚度为50 ～ 250 m，其优质的天然条件为该地区地下水的开发利用提供了保障（图2（a））。典型区2019年3月地下水流场如图2（b）所示。

图2 三江平原典型区富水性分区及2019年3月地下水流场

1.2 数据来源 选取的遥感影像数据［16］源自中国科学院资源环境科学数据中心（https：//www.resdc.cn），包含1980年的Landsat-MSS遥感影像数据，1990年、2000年、2010年的Landsat-TM/ETM 遥感影像数据以及2018年的Landsat 8 遥感影像数据。根据徐新良等［17］的中国土地利用/覆盖数据分类体系，考虑到三江平原典型区实际地物分布与遥感影像的可解译性，将该范围内的土地利用类型分为：水田、旱地、水域、林地与草地、城镇用地、农村居民点、其他建设用地及未利用土地［18］（表1）。

表1 土地利用类型分类

利用ArcGIS 10.6软件进行遥感影像数据处理，并绘制1980、1990、2000、2010及2018年五期的土地利用类型分布图（图3）。

图3 三江平原典型区1980—2018年五期土地利用类型分布

2 研究方法

对土地利用结构与时空变化规律的分析主要体现在土地利用类型的数量变化和动态变化等方面，土地利用结构的稳定直接影响流域水均衡状态的稳定。

2.1 土地利用变化研究

2.1.1 土地利用动态度与分形维数

（1）土地利用动态度。土地利用动态度是指研究区在某段时间内某土地利用类型和数量的变化情况，可以反映区域土地利用类型的变化剧烈程度［19］。

式中：K为某种土地利用类型的动态度；St1、St2分别为研究初期与末期某种土地利用类型的面积，m2；T为研究时间间隔，a。

（2）分形维数。分形维数［20-21］是对某一结构的自相似性特征进行量化表达的指标，能够反映一个分形体的不规则程度。利用分形维数既可以对该结构的复杂性进行定量测度，也可以对其稳定性进行表征参量。研究土地利用空间结构的分形特征，主要是定量描述其面积的大小及其边界线的曲折程度，基于已知的各土地利用类型的周长、面积等信息，采用周长—面积法拟合出面积关于周长的双对数回归分析模型，其计算公式为：

式中：A(r)为斑块面积，m2；D为分形维数，一般情况下，D∈[1，2]，D越大，该类用地规则程度越低，其结构越复杂；P(r)为斑块周长，m；C为常数。

基于分形维数D计算得到稳定性指数S，并进行显著性检验，确定系数R2均大于0.9 则结果可信，其计算公式为：

式中：S为稳定性指数，S∈[0，0.5]，S越大，该类用地越趋于稳定状态；D为分形维数。

2.1.2 标准差椭圆法 标准差椭圆法［22-24］（SDE）是进行数据点集或区域的中心趋势、离散和方向趋势分析的最直接方法，该方法通过以中心、长轴、短轴、方位角为基本参数的空间分布椭圆定量描述研究对象的空间分布整体特征。重心计算公式为：

2.2 大气降水入渗计算分区 根据入渗补给区岩性的不同对研究区进行分区，选取不同的大气降水入渗补给系数（图4及表2），从而计算得到大气降水入渗补给量。

表2 大气降水入渗补给系数分区

图4 大气降水入渗补给系数分区

3 结果与分析

3.1 近40年土地利用变化特征分析

3.1.1 土地利用动态分析 三江平原典型区1980—2018年五期的各土地利用类型变化情况（表3）显示，水田、城镇用地与农村居民点变化面积在各研究时段内均为正值，表明其面积在近40年内呈逐年增加的趋势，其中水田面积在1980—2000年间增幅稳定，2000—2010年增幅最大，面积增加了4377.72 km2，动态度达到峰值58.97%，2010—2018年水田面积增加减缓，动态度迅速下降；城镇用地面积在1980—1990年增幅最大，1990—2000年面积增加最为缓慢，动态度仅为0.22%，2000年之后，面积增幅虽大于1990—2000年，但增速减缓，动态度直线下降；农村居民点面积在1980—1990年增幅最大，动态度达到峰值8.47%，从1990年起，面积增加趋势减缓，动态度在0.01%～0.61%之间波动，面积变化不明显。

表3 三江平原典型区1980—2018年各土地利用类型变化情况

旱地面积近40年呈现先增加后减少的变化趋势，1980—2000年动态度直线下降，2000—2018年动态度趋于平缓，部分旱地逐渐变为水田；林地、草地与未利用土地的变化面积在各研究时段内均为负值，表明林地、草地与未利用土地的面积呈逐年下降的趋势，且均在2010—2018年下降幅度最大，动态度分别为-4.71%与-6.30%，结合图3来看，大面积的林地、草地与未建设用地先变为旱地再变为水田。水域面积在1990—2018年间呈逐年减少的趋势，降幅逐年增加。其他建设用地面积自1990年开始变化，呈现先减少后大幅增加的趋势，2010—2018年动态度达到峰值148.02%。

由于影响地下水位埋深及流域水均衡要素的主要为水田和城镇用地，故本文主要对这两种土地利用类型进行分析与讨论。对三江平原典型区1980—2018年五期的水田和城镇用地斑块数、斑块周长、斑块面积进行统计（表4），并得到各土地类型的双对数回归分析模型与分形维数（表5）。由于近40年水田与城镇用地的R2均达到0.9以上，表明这两种土地类型的面积与周长存在显著的线性相关关系，分形维数具有一定意义。

由表4、表5及图3可知：近40年来，水田与城镇用地都存在较为显著的分形特征。两种土地利用类型的分形维数均先增大后减少，表明它们的结构先逐渐复杂，之后其结构异质性和景观破碎程度不断下降，各边界形状逐渐趋于平整，结构及变化趋势逐步趋于稳定。

表4 三江平原典型区1980—2018年部分土地利用状况

表5 三江平原典型区1980—2018年部分土地利用类型分形维数及稳定性指数

水田是研究区近40年来面积变化最大的土地利用类型，自1980年起在研究区面积占比逐渐增大且扩展面积最大，1980年水田面积仅占总土地面积的0.47%，而到2018年时，水田成为研究区内第一大用地类型，其面积占总土地面积的47.68%。从近40年分形维数与稳定性指数来看，水田的分形维数从1980年的1.0965 升高至2000年的1.2367，2000年后开始下降，到2018年分形维数降至1.1527，稳定性指数呈先降低后升高的趋势，表明水田在研究初期即1980年由于面积较小，复杂性最低，稳定性最强，随着工农业的迅速发展与城市化进程的加快，水田面积不断扩大，城镇及交通用地空间不断扩张等，导致城镇周围的水田不断转变为建设用地，水田的空间形态也受到了一定程度的破坏，到2000年时水田的复杂性达到最强，稳定性处于最低。随着一系列的农业生态保护工程与措施的推进以及城市化的进一步发展，水田空间规则性与稳定性得以恢复和提高，到2018年水田空间结构的稳定性已维持在一个较好的状态中，农业发展稳定。

城镇用地面积变化较小，虽在近40年中面积逐渐增大，但在总土地面积中占比较小，维持在0.05%～0.39%之间。1980年城镇用地斑块数多，但平均斑块面积小，表明城镇用地形态细碎，密集散布在研究区内，主要分布于西南部。1980～1990年城镇用地斑块数急剧减少，但1990年城镇用地面积较1980年增加了18.69 km2，表明城镇用地的扩张在研究时段开始的十年内十分活跃，该时期城市化进程加快，零星城镇用地得以合并，分形维数达到研究时段内最大值（D=1.5321），稳定性指数仅0.0321，土地复杂性最高，城镇用地自1980年起扩展后，由于还在发展阶段，空间形态不规则，空间布局不合理。2010年之后，城镇用地斑块数稳定在23个，分形维数稳定在1.15左右，稳定性指数也稳定在0.34 左右，但面积依然持续扩张，表明城镇用地布局趋于合理化，扩建后形态趋于规则，当前城市化发展已趋于稳定，围绕现有的城市中心稳步发展。

3.1.2 土地利用时空迁移分析 对1980—2018年五期土地利用类型影像进行空间叠置分析，分别将水田和城镇用地的动态演变阶段的重心进行叠合，形成两种土地利用类型的重心迁移图（图5），并计算得到其重心迁移距离（表6）。由图5 与表6 可知，研究时段内水田和城镇用地重心的迁移方向主要为西南—东北—西南，两种土地利用类型的重心均在1980—1990年时迁移距离最大，分别为20.41 km和56.06 km。

表6 1980—2018年土地利用类型重心迁移距离及迁移方向

图5 1980—2018年水田及城镇用地重心迁移图

水田重心迁移轨迹可分为四个阶段：第一阶段为1980—1990年，重心由西南向东北迁移，该时段东西方向迁移幅度强于南北方向，迁移速度最快；第二阶段为1990—2000年，重心继续向东北方向迁移，北移趋势明显但迁移速度缓慢，仅迁移了6.33 km；第三阶段为2000—2010年，重心向东迁移；第四阶段为2010—2018年，重心由东北向西南方向回移了19.95 km，迁移速度加快，向西回移趋势明显，到2018年时重心落在了1980年与2010年的中间位置。整体来看，研究时段内水田重心迁移幅度较小，均在整个研究区的中心地带（七星农场五十五队—前进农场十四队）移动。

城镇用地迁移轨迹可分为三个阶段：第一阶段为1980—2000年，重心由西南向东北迁移，该时段迁移距离为83.64 km，其中1980—1990年以5.61 km/a 的速度加速迁移，1990—2000年间迁移速度降至2.76 km/a；第二阶段为2000—2010年，重心以4.11 km/a 的速度加快向西南方向回移了41.14 km；第三阶段为2010—2018年，重心仅向南偏东方向迁移了0.07 km，迁移距离最小，重心位置自2010年起基本稳定，城市化进程以该处为中心向周围推进。整体来看，近40年城镇用地重心迁移较水田重心迁移跨度大，涉及范围广。

3.2 LUCC对地下水位埋深的影响分析 由于所收集资料有限，本文只讨论2000—2018年三江平原典型区土地利用类型变化对地下水位埋深的影响。由图6与表7可知，2000—2018年间研究区地下水位埋深呈由东北向西南方向逐渐增大的趋势。结合图3与图7可知，2000—2018年研究区降水量增加了49.60～259.19 mm，虽然在富锦市与同江市附近降水量增幅最大，但富锦市与同江市附近的地下水位埋深并未相应地减小，这是因为地下水位埋深影响较大的因素是地下水开采（对应地区为水田分布区），研究区东北部在2000—2018年间林地与草地逐渐转化为旱地，且距离河流较近，地表水资源丰富，部分地区采用地表水进行农田灌溉，对地下水的需求较小，故地下水位埋深减小，抚远市附近地下水位埋深减小最为明显，地下水位埋深减少了3.86 m。

表7 2000—2018年不同地下水位埋深变化范围的各土地利用类型面积变化量

图6 2000—2018年地下水位埋深变化分布

图7 2000—2018年降水量变化分布

地下水位埋深增加小于1 m 的区域面积为13 341.44 km2，占整个研究区面积的63.85%，且主要集中于研究区中部，西南部也有小部分面积涉及。到2018年时，研究区中部分布有大面积水田，且65.35%的城镇用地和部分工业园区也集中在该区域。人口密集，用水量大，是该区域地下水位埋深增大的主要原因之一。2000—2018年该区域水田面积增加了6891.95 km2，城镇用地面积增加了37.74 km2，旱地面积减少了3840.31 km2，林地与草地面积减少了1478.13 km2，未利用土地面积减少了1355.01 km2。土地利用类型面积的增减量达到了地下水位埋深下降小于1 m 的区域总面积的1/2，表明该区域大面积扩张水田和城市化的快速发展对地下水资源进行了一定量的开采。

地下水位埋深增加超过1 m 的区域面积为4745.40km2，占整个研究区面积的22.71%。该区域在2000—2018年间水田面积增加了1877.92 km2，城镇用地面积增加了7.07 km2，未利用土地面积减少了1090.09 km2。该区域地势较高，地下水经开采后补给恢复较慢，因此，该区域的地下水位埋深在整个研究区最大。

研究区自西南向东北地势逐渐降低，沿东北方向地下水补给恢复加快，东北部多条河流汇集，灌溉对地下水需求减少，再伴有区内降水分布的不均匀，使得地下水位埋深自西南向东北逐渐减小。

3.3 LUCC对水均衡状态的影响分析

3.3.1 城市化进程对大气降水入渗补给量的影响 1980—2018年大气降水入渗补给量变化情况见表8。由表8可知，由于研究时段内城市化进程加快，不透水面积增加，部分降水成为地表径流，降水入渗补给面积减少，从而影响大气降水入渗补给量。1980—2018年不透水面积的增长速度先降低后增加，2000—2010年不透水面积以3.35 km2/a的速度迅速增加，2010—2018年次之。虽然2018年的不透水面积为1980—2018年间最大值，但由于该年降水量较大，大气降水入渗补给量并不为近40年来的最小值。这是因为城镇用地面积在整个研究区占比较小，不透水面积变化程度对大气降水入渗补给量影响也不显著，导致研究区大气降水入渗补给量的大小更取决于当年降水量的多少。随着城市化进程的推进，若城市继续扩张，将会对研究区水均衡状态产生更大的影响。

表8 1980—2018年大气降水入渗补给量变化情况

3.2.2 耕地面积变化对灌溉回归补给量的影响 根据所收集的资料，本文仅讨论2010年与2018年的耕地（水田与旱地）面积变化对研究区灌溉回归补给量的影响。由表9与图3可知，2010—2018年耕地面积增加3828.17 km2，其中水田面积增加5383.79 km2，旱地转化为水田的面积为3720.98 km2，灌溉回归补给量增加2.40亿m3。可见，虽旱地2018年的面积较2010年减少了-1555.62 km2，但水田面积的大幅增加，使得灌溉回归补给量大幅增加，故在灌溉回归补给量的计算中，水田面积的变化是其主要影响因素。

表9 2010—2018年灌溉回归补给量变化情况

4 结论

本文基于三江平原典型区1980—2018年间的土地利用类型数据，利用动态度、分形理论与标准差椭圆法分析了研究区该时间段内土地利用变化动态特征与时空迁移特征，并在此基础上探究了土地利用变化对地下水位埋深与流域大气降水入渗补给量和灌溉回归补给量的影响，得出以下结论。

（1）1980—2018年间研究区土地利用类型发生了较大的变化，其中水田、城镇用地与农村居民点面积呈逐年增加的变化趋势，旱地面积与其他建设用地面积呈先减少后增加的变化趋势，林地、草地与未利用土地面积呈逐年下降的变化趋势，水域面积呈先增加后减少的变化趋势。

（2）在土地利用动态方面，1980—2018年间来水田与城镇用地的面积和周长存在显著线性相关关系与分形特征。两种土地利用类型分形维数分布在1.0～1.6之间，均先增大后减少，稳定性指数由强—弱—强。水田是研究区近40年来面积变化最大的土地利用类型，在2018年成为研究区内第一大用地类型。城镇用地面积变化较小，2010年后布局趋于合理化，扩建后形态趋于规则，城市化发展趋于稳定。

（3）在土地利用时空迁移方面，近40年来水田与城镇用地重心发生迁移，且迁移方向主要为西南—东北—西南，水田重心迁移幅度较小，均在整个研究区的中心移动；城镇用地重心迁移跨度较大，重心自2010年起基本稳定，城市化进程以该中心向周围推进。

（4）土地利用类型的变化是影响研究区地下水位埋深与水均衡状态的主要原因之一。研究区地下水位埋深由东北向西南方向逐渐增大，中部地区分布的大面积水田对地下水需求较大且该区域地势较东北高，地下水补给恢复较慢，使地下水位埋深有不同程度的增大。此外，城市化发展对大气降水入渗补给影响较小，但若城市不断扩张，将会对大气降水入渗补给产生较大的影响；水田面积的变化是灌溉回归补给量变化的主要影响因素，2010—2018年水田面积增加5383.79 km2，灌溉回归补给量增加2.40亿m3。