利用gis预测和降低地面沉降的时空发展

采用GIS方法建立地下水流和地面沉降的时空预测系统，并应用于沉降持续发展的日本沿海平原。将与地下水流和地面沉降相关的不同基础数据进行数字化，并输入到GIS数据库。采用GIS对整个平原的地表水文地质循环进行模拟，据此获得了地下水入渗量。通过GIS数据库将数据转化为地下水流动代码(MODFLOW)，建立三维地下水流动模型，完成过去21年的非稳定地下水流动模拟。最后，用VB编写了含水层和隔水层的地面沉降的程序。根据不同的地下水抽水量对未来的地面沉降进行预测，研究结果对于制定减缓地面沉降的决策提供了重要信息。一、概述在世界的许多地区都出现了地面沉降，特别是在人口稠密的三角洲地区，造成了巨大损失。目前，世界上有150个以上的地区，包括墨西哥、中国、台湾、意大利、日本和美国(德克萨斯、加利福尼亚、亚利桑那、路易斯安那、内华达州、爱达荷州和其它州)，都出现了地面沉降。最严重的地面沉降是由人为因素诱发的，如过度开采地下水、石油、天然气和地热流体。地下水开采造成的地面沉降幅度在几毫米到10米以上，沉降范围与地区有关，一般从几平方公里到10000平方公里以上，对自然和社会环境造成了直接和间接影响。根据统计资料（Johnson，1991），估计全世界由于地面沉降造成的破坏和修复费用为数十亿美元。如在美国休斯顿，每年的地面沉降费用为3170万美元（Jensen，1985）。在最近几年，地面沉降已经成为一些发展中国家严重的地质环境问题，在这些国家，对自然资源（特别是地下水资源）的需求在增加，经济发展和环境保护的矛盾也越来越明显。典型的情况是在中国，由于大规模的抽取地下水来维持迅速发展的经济，约有50个城市出现了地面沉降。过去采用了许多方法来解决这一人为诱发的环境地质问题，以前的地面沉降调查分为两类：1．在特定点进行调查，主要是根据地质钻孔资料、水位变化和地面沉降量，不考虑地面沉降的空间分布，来获取详细的信息（Helm，1975，1994；Leak，1991等）。2．进行区域地面沉降调查，模拟区域地下水流动和地面沉降，并提出控制地面沉降的地下水抽水方案（Esaki等，1995，1996；Teatini，1995）。根据过去的调查，地面沉降是非常复杂的系统问题，包括不同的空间分布信息，如地质、地形、土地利用、降水、蒸发、含水层系统、水文地质参数分布、抽水、地下水流、入渗等。要对这样一个系统进行模拟，要采用多种方法和系统的调查研究。然而在过去的大多数情况下，对这一系统各组分进行的是孤立的研究，例如，土地利用调查一般为制定土地利用计划，对蒸发蒸腾的研究的目的是为农业发展。在将所有组分结合起来分析时，量化过程和分析大量的空间分布数据是非常重要的，如分析水力传导率、地形、入渗的空间分布等。传统上，所有这些数据处理都是手动的，这样工作效率很低，而且建立模型的可靠性和分析的精确性都有疑问。另外，由于地面沉降一般都是在长时间连续发生的，也要考虑时间特性和环境对它们的影响。在最近一些年，具有多功能强大的处理、管理和分析空间分布数据的地理信息系统（GIS），在许多领域，包括在水资源管理和自然灾害缓解方面，引起了人们的高度重视。GIS系统可以进行数据采集、输入、处理、转化、目视、合并、查询、分析和模型化。GIS在进行空间分布数据的处理和分析过程中具有强大的功能，采用GIS也可以有效地进行自然灾害的环境影响评价。在GIS基础上用多种方法将地面沉降进行模型化。GIS作为处理与地面沉降有关的空间数据的工具和平台，将分散的资料，如地表水水文地质循环、地下水流动和地面沉降模拟结合起来进行分析。通过日本沿岸平原的特定事例研究，介绍建立空间和时间预测系统的GIS程序。二、日本和Saga平原的地面沉降长期以来，日本的地面沉降就是一个环境问题。自上世纪30年代以来，日本城市和工业迅速发展，为满足工农业发展和日常生活用水等，对地下水资源的需求持续增加。在许多地区过度开采地下水，地面沉降成了一个全国性的环境地质问题。自1979年，在大工业城市，地面沉降的总趋势在很大程度上减慢了，这是由于政府机构和工业水法限制了地下水开采。但是，在近几年，与大城市地面沉降减缓相反，一些中等城市和平原（如Saga平原，Nobi平原等），仍受地面沉降的影响，主要是为满足农业用水进行地下水开采。Saga平原位于日本西南，是一典型的沿海低地，平均地表海拔2.41m。在涨潮时几乎整个平原都在水位下。自1957年来平原一直发生地面沉降，目前，最大累积地面沉降量为124cm，地面沉降区范围为320km2，其中沉降量大于10mm的面积为191km2（1999年资料）。三、Saga平原的水文地质条件和地下水抽水情况（一）水文地质条件Saga平原北部为MesozoicSefuri山脉，东部为由Sangun变质岩组成的Chikuhi山脉，西部是Kishima山，南部是Ariake海。河系经过平原最终汇入大海。平原分为三个区域：西部的Shiroishi区，中部的Saga城区和东部的Chikugo区。Saga平原下分布有厚度为200m的沉降物。这些沉降物分为A、B、C、D、E和F地层（Oshima，1988）。地层A是Ariake冲积粘土沉积物，主要由软粘土和粉土组成，局部夹砂层。Saga平原的Ariake粘土层厚度在10～30m，平均厚度为20m。地层B是洪积海相沉积物，主要由砂组成。地层C是浮石/泡沫状火山灰。地层D、E和F都是由海相砂和粉土组成。地层B、C、D、E和F是平原的含水层，主要是从这些含水层抽取地下水。在B、C、D、E和F地层中，分布有许多粘土层，这样Saga平原的含水层系统就更为复杂。（二）地下水抽水、水位波动和地面沉降最近几年，Saga城区的年抽水量为300万m3，在Shiroishi地区（Saga平原西部），抽水量是700万m3。Saga城区和Shiroishi地区的抽水原因完全不同。根据1994年的统计资料，Saga城区84%的抽水量用于工业发展，而在Shiroishi地区80%的抽水量用于稻田灌溉，而且在夏季稻田灌溉用水量增加，地下水抽水量也较大。根据1975～1997年在Shiroishi地区观测的水位波动和地表高度变化资料，可以看出地下水水位变化与地下水抽水量有关，在夏季地下水位下降，在其它季节水位上升。同时，地表高程也是在夏季下降，在其它季节由于残余压缩有所回弹。据此可以得出为模拟地面沉降，有必要对整个平原区，包括海水相对较浅，仅为2～3m的Ariake海和海底50km厚的Ariake粘土层的地下水流动情况进行模拟。四、GIS方法如前所述，对整个平原区的地下水流动机制进行分类，包括补给、排泄、径流和水位变化等，是模拟地面沉降的重要步骤。Saga平原的地下水系统由非承压水和承压水组成，非承压水主要是存储在Ariake粘土层的浅层，承压水存储在B、C、D、E和F洪积层中。在承压含水层中进行季节性抽水，水位变化与抽水量有关。承压水主要由平原边界的水平流补给，砾石组成使得地表水易于入渗，在垂直方向连续分布有20m厚的Ariake粘土层，在大多数地区成为隔水层。在大多数地区，浅层潜水水位在地表下0.5m（日本环保局，1982），水位变化受地下水入渗的影响。在应用三维非承压－承压地下水流动模拟时，有必要进行水文地质模拟来判断地表水入渗情况。在进行地表水水文地质循环模拟和进行地下水流动与地面沉降模型化时，要对大量的空间数据进行存储和处理，采用GIS作为处理和分析空间数据的平台和工具，GIS进行地面沉降的空间和时间预测系统见图1，由以下三部分组成：\o"点击图片看全图"\o"点击图片看全图"五、地表水循环模拟（一）模拟方法一般而言，有两种地表水循环模拟方法：Lamped模型方法和分布模型方法。Lamped模型方法包括水库模型（或水箱模型）和径流系数模型。分布模型能够考虑到地表、坡角和植被等的影响，但是建立分布模型很难而且耗时，本研究采用了径流系数模型。径流系数方法要考虑的最基本问题，是地表水接受的降水量（作为输入）和从河流，地下水和由于蒸发蒸腾、径流和向地下入渗等造成的自身消耗的农业灌溉用水（作为输出）。输入和输出应当达到平衡，地表水水文地质循环模拟公式如下：（1）其中，Ra是降水量，Ev是蒸发蒸腾量，Ru是径流，Ar是农业灌溉用水，In是入渗量。根据上式，在其它条件已知的情况下可以得到入渗量。收集与降水、农业灌溉、蒸发蒸腾和径流相关的资料，将其空间分布作为GIS数据层。采用GIS数据处理（层处理），可以获得入渗的时空分布。用ArcGIS8.1软件，将所有GIS数据格式转化为矢量数据或格栅数据格式。（二）降水在Saga平原，有8个国家级气象观测站，观测项目包括降水量、温度、风速和其它，收集了1979～1999年的资料。将观测站的位置以点的形式输入到GIS，通过GIS的嵌接处理，将降水数据以Excel的形式与观测点的属性数据表连接。在GIS环境下点击一个观测站点，就可以迅速确定过去21年（1979～1999年）任意一天的天气情况（降水、温度等）。根据这些观测点，采用GIS宏语言编Thissen程序，获得每一个观测点的Thissen多边形，而且观测点具有相同的属性表。在同一Thissen多边形区，认为降水量是相同的。在Saga平原，由于地下水抽水量统计是根据1km×1km的网格进行的，因此地下水流动和地面沉降模拟也应当根据这样的网格进行。因此，对于地表水水文地质循环模拟，最终的目的是计算每1×1km的网格中的降水量。通过将1km×1km的网格多边形和Thissen多边形叠加，Thissen多边形的属性数据可以转化为网格组成。对于那些与Thissen分割线相交的网格组成，通过对相交网格组成进行面积加权确定天气状况。根据以上的GIS处理，以1km×1km的网格形式，获得1979～1999年的降水空间分布。（三）蒸发蒸腾通常采用Thornthaite方法进行蒸发蒸腾量计算，Thornthaite方法计算公式如下：（2）其中，Et是可能的蒸发蒸腾量（单位cm/月），T是月平均温度，，。当已知一点的月平均温度时，可以通过以上等式计算蒸发蒸腾量。根据以上分析，气象观测站是GIS点数据形式，温度记录是属性数据形式。采用温度数据处理GIS属性数据表，根据以上公式可以计算每一观测点的月蒸发量。采用相同的程序（建立Thissen多边形，并与网格多边形叠加），以1km×1km网格多边形形式，获得1979～1999年蒸发蒸腾量的空间分布。（四）径流在降雨时，一些雨水流经地表最终汇入河流，其余的雨水流入地下，通常称为径流量（Ru），可以根据下式进行计算：（3）其中Ra是降水量，是径流系数。根据式（3）计算径流量时，必须知道土地利用情况的空间分布，再对不同的土地利用类型确定径流系数。日本地质调查局提供了1971、1991和1999年每年全国土地利用数字图。这些图由将土地利用类型作为属性数据的100×100m多边形组成。在本项研究中用到了1991年的数字地图。根据2000年日本建委出版的“侵蚀控制技术标准”，确定不同土地利用类型的径流系数，见表1。通过将土地利用图属性表和径流系数表进行空间嵌接处理，将径流系数加入到土地利用多边形属性表中。根据获得的降水分布图层、土地利用图层和式（3），用GIS按月份计算1979～1999年每100m×100m多边形的径流量，计算出的径流量作为土地利用多边形中新的属性。将100m×100m多边形叠加在1km×1km网格多边形上，以1km×1km网格多边形形式，获得1979～1999年不同月份径流量的空间分布。（五）农业灌溉用水在Saga平原，大部分地下水和一些河水用作稻田的灌溉用水，稻田占整个平原地区用地的70%。一般来说，农业灌溉用水（Ar）可以根据下式得出：（4）其中Dh是在灌溉季节的稻田耗水量，A是灌溉面积。根据过去在Saga平原许多稻田的野外观测，Dh平均值为15mm/天。根据数字土地利用图，采用GIS提取法得出稻田分布图。根据稻田分布图和式（4），计算每100m×100m稻田多边形的农业灌溉用水量，并得出灌溉季节（6～9月）的灌溉用水空间分布。采用相同的程序，重新计算1km×1km网格多边形形式的农业灌溉用水量。表1不同土地利用类型的径流系数土地利用径流系数土地利用径流系数稻田0.70交通用地0.85小块地0.52水库1.00果园0.63河流0.75森林0.52海滨0.20弃耕地0.60海域1.00建筑用地0.90--（六）入渗量GIS处理根据以上的GIS空间数据处理，获得4个GIS层：降水、蒸发蒸腾、径流和农业灌溉用水。所有的图层都采用相同的1km×1km多边形，共有973个多边形。每一层有一个属性表，有252行，按月存储1979年1月至1999年12月的相关数据。为获得入渗量空间分布，用到了4个属性表中的数据，见下式：（5）其中In是入渗量。根据上式可以得出1979～1999年每一网格多边形的每月入渗量。根据计算结果，可以看出在灌溉季节，平均入渗量约为10mm/天，在其它季节不足1mm/天，这一计算结果与Tohara和Kato于1978年进行的野外入渗观察结果一致。六、三维地下水流模拟（一）水文地质概念模型如前所述，从地表下20m有一连续分布的Ariake粘土层，成为Saga平原的隔水层。在Ariake粘土层下有一厚度约为200m的洪积层，地下水主要是取自该地层。水文地质模型见图2。根据下式进行三维地下水流动模拟：（6）其中h是水头（m），k是渗透率（m/天），Ss是单位产水量（l/m），Qd是抽水量（m3/天），Qr是垂直补给（m3/天）。用MODFLOW模拟三维非稳定地下水流。中，在所有点之间进行内插，这样可以以光栅数据的形式获得隔水层的空间分布。根据隔水层厚度和地表高程光栅数据，经过光栅数据处理可以获得隔水层底板的空间分布。采用与地形数据处理相同的数据输出程序，将隔水层底板的空间分布转化为MODFLOW代码。（五）抽水量空间分布Saga地区连续多年对地下水抽水量进行了统计