0 引言

羊庄盆地岩溶水系统是鲁中南山区较为重要的岩溶水系统，补给来源充沛，岩溶水丰富，调节功能优良。目前羊庄盆地岩溶水的开发利用已达到较高水平，区内先后建立了后石湾、西石楼、魏庄、羊庄等大中型水源地。水源地的开发为当地经济发展及人民生活水平的提高做出了重要贡献，但也对生态环境造成了不良影响。目前，对岩溶地下水资源进行有效管理已成为热点课题，而地下水管理问题通常通过模拟-优化方法来解决^[1]。模拟模型能解决在一个给定的开采方案前提下的预测问题，但它不能确定哪个开采方案能获得最优结果; 而通过优化模型则可确定出满足设定目标函数和约束条件下的最优开采方案。因此，模拟和优化是目前对实际地下水系统进行机理分析，模拟预测和开采方案优选的最佳方法^[2-4]，其最大的优点就是既可以模拟复杂的地下水流系统，同时还可得到满足目标函数和约束条件下的最优开采方案。目前国内外大多采用REMAX、MODOFC、GWM、MODMAN、MGO等地下水管理软件，将水流模拟模型和优化模型耦合，从而确定出地下水资源的优化开采方案^[5-11]。随着地下水管理模型软件的逐渐成熟和发展，耦合方法已不再是影响优化计算的最大障碍，而数值模拟的可靠性则成为影响优化结果的关键因素，只有在模拟模型准确的前提下，才能保证最优化的过程遵循实际地下水系统的运动规律(地下水运动规律由模拟模型来呈现)^[12]。而模拟模型中源汇项的准确与否直接关系着数值模拟的可靠程度，源汇项中最难确定的为分散的农业灌溉和预测阶段的降水量。长期以来模型计算者大多利用灌溉定额法估算农业开采量，选取不同保证率下的降水量作为未来降雨量，过大的源汇误差严重地影响了数值模拟的可靠程度，同时也影响着优化结果的准确性。

本文以山东羊庄盆地为例，在利用“作物需水量”方法对农业开采量进行精确化计算的基础上，首先根据已有降水数据基于非平稳时间序列对未来的降水量进行准确预测，摆脱了定额法确定农业开采量和人为构造降雨序列的粗糙性和随意性；然后将校正、验证后的地下水流数值模型与优化模型相耦合，同时借助通用的优化模型软件Ground-Water Management Process(GWM)，建立了羊庄盆地地下水资源优化管理模型；最后以开采量最大化为目标，以水位不低于限定的控制水位为约束条件，对盆地开采潜力进行了评价; 进而提出各个集中水源地最优取水方案，以期能为羊庄盆地岩溶地下水资源合理开发和管理规划提供科学依据。

1 研究区水文地质条件

羊庄盆地位于山东省南部，包括滕州市的羊庄镇和枣庄市山亭区的部分地区。在地貌上属于鲁中南低山丘陵区，为一四周环山的向斜盆地，面积约650 km²，标高50~600 m。

羊庄盆地岩溶水系统是鲁中南山区较为重要的岩溶水系统，其四周较封闭、出口狭窄、岩溶发育，具有完善的输导、调蓄功能，存在基本统一的水动力场，是典型的封闭式泉排型岩溶水系统。盆地内部广泛分布着碳酸盐岩类岩层，岩层均从四周向腹地倾斜。其西北、北及东部由寒武系灰岩、古老变质岩和火成岩组成地表分水岭，南及西南部通过枣庄断裂和化石沟断裂与煤系地层接触，是一个较为封闭的岩溶蓄水盆地。新薛河是区内唯一的地表河流，河谷深切，多处河段与地下水存在互补互排关系，是地下水获得补给与溢出排泄的主要通道之一。羊庄—西集一带的盆地腹地第四系广布，并发育一条长15.4 km，宽375~1 500 m的第四纪古河道，砂砾石层直接覆盖于下伏灰岩含水层之上，岩溶水可通过顶托补给孔隙水的形式向区外径流排泄。

羊庄盆地由长龙断裂和曹王墓断裂将盆地分割成辛召断块、山亭断块和羊庄断块3个水文地质子单元(图 1)。辛召断块和山亭断块两个子单元位于羊庄盆地的中、上游地区，地势较高，地形起伏较大，地下水的循环径流条件和含水介质的富水性均相对较差，地下水资源较为贫乏；而羊庄断块子单元处于羊庄盆地下游，地下水的补给、循环条件优越，含水介质的富水性强，是羊庄盆地岩溶地下水循环径流和富集储存的主要场所(图 2)。

2 水文地质概念模型 2.1 地下水系统边界条件的概化

羊庄盆地为相对独立、有较完整补、径、排条件的向斜盆地式封闭岩溶含水系统。盆地西北及东部边界由不透水的泰山群变质岩、中下寒武系页岩夹灰岩等岩性组成；南部通过陶枣断裂与煤系地层相接，使其两条断裂本身即构成隔水边界；只有三山头—龙山头之间有一狭窄的地段，其地势低、地形标高小，成为整个盆地地表水、地下水的唯一出口^[13]。

1) 侧向边界的概化：为避免人为边界对模型计算造成的困难，在充分利用盆地自然边界的基础上对全盆地地下水流系统进行整体模拟。模拟区平面上，除整个羊庄盆地出水口处(三山头—龙山头断面)属第二类定流量排泄边界外，其余均为隔水边界(图 1)。

2) 底面边界：据区内钻孔资料，含水层大多集中在埋深220 m以上，以下岩溶基本不发育，故将埋深220 m位置概化为隔水边界。

2.2 源汇项的概化

模型的主要源汇项包括：新薛河及其支流的补排量、农灌回归入渗补给、大气降水入渗补给、泉水排泄和地下水的人工开采量。本文着重讨论对农业分散开采进行精确定量化的新方法，其他源汇项水量的确定按照一般常规方法进行。

对于分散的农业灌溉，由于灌溉井离散分布的特点，很多水井不受灌溉管理部门的直接管理，农户在地下水的使用上具有很大随意性，因此，管理部门所提供的农业灌溉开采量统计值往往偏小^[14]。笔者此次利用灌溉面积和作物净灌溉需水量来估算地下水开采量。确定评价区农业分散开采量主要分为作物净灌溉需水量的确定和评价区耕地灌溉面积的估算两步。

1) 作物需水量的确定

作物需水量(T_Ec)可由参考作物蒸腾发量(T_E0)和作物系数(K_c)来确定，其计算公式^[15]为

T_E0可采用联合国粮农组织推荐的彭曼公式^[15-17]计算：

式中：Δ为饱和水汽压与气温曲线的斜率，kPa/℃；R_n为冠层表面净辐射，MJ/(m²·d)；G为土壤热通量，MJ/(m²·d)；γ为湿度计常数，kPa/℃；T为平均气温，℃；u₂为2 m高处的风速，m/s；e_s为饱和水汽压，kPa；e_a为实际水汽压，kPa。

作物系数在一定范围内为常数，因此可参考邻近站点K_c值，此次参考鲁西南地区莱河代表站的取值，见表 1。羊庄盆地农作物以冬小麦和夏玉米为主，兼有少量蔬菜和果树，为了计算作物需水量时方便选取合理的作物系数，此次评价假设整个盆地农作物只有小麦和玉米两种。

根据公式(1)和(2)，参考FAO－56中^{[10, 12]}参数的详细计算过程，采用羊庄气象站气象资料，计算出评价区两种作物2013—2016年逐月作物需水量，见图 3。

计算时忽略降水过程中的损失，假设有效降水量等于降水量，用P代表降水量(mm)、T_Eg代表净灌溉需水量(mm)，则当T_Ec＞P时，T_Eg=T_Ec-P，当T_Ec < P时，T_Eg=0；最终确定出评价区两种作物2013－2016年逐月作物净灌溉需水量，见图 3。

2) 灌溉区面积的估算

基于野外实地考察所建立的遥感解译标志库，笔者采用人工解译与机助解译相结合的方法^[19]，对整个羊庄盆地耕地面积进行了遥感解译，通过解译得到羊庄盆地耕地面积为55.2 km²，空间分布见图 4。

3) 农业灌溉开采量估算

评价区农业灌溉总需水量由作物净灌溉需水量和遥感解译的耕地面积相乘，确定为平均10.3万m³/d；而农业灌溉开采量应将渠灌水量从农业灌溉总需水量中减去。

对于渠道引水灌溉，因为灌溉部门有统一规划的灌溉方案，因此统计数据较为真实。据水管部门提供的渠道灌水统计资料可知，研究区内渠灌引水量为5.43万m³/d，则农业灌溉开采量平均约为4.86万m³/d。

鉴于盆地内农灌井分布太多，在模型中无法逐一设置，因此我们将农业灌溉开采量按照农田面积平均分配到农田所在单元格中，作为面状(负补给)开采加入模型中^[20]。

3 模型识别及结果

此次模拟选取2013-05-26—2016-09-26作为模型的识别验证期，划分为41个应力期，其中2013-05-26—2015-01-06的资料用于模型校正，2015-01-07—2016-09-26的资料用于模型检验。以观测孔的地下水位动态作为拟合对象来检验模型^[21-22]。

首先根据区域水文地质条件、水文地质勘探结果等初步划分各参数分区，给出各参数初始估计值和变化范围, 经过多次调试、优选, 最终使各观测孔的模拟水头值与实测值达到最佳拟合^[23-24]，确定将全区共分为19个含水层参数分区以及8个降雨入渗补给系数分区(图 5，6)。各分区参数值见表 2和表 3。总体来看，模拟得出的水文地质参数与抽水试验求得的参数具有一致性，参数取值较为合理。

由拟合结果(图 7)可看出，研究区25个水位观测孔拟合效果较好，且各测点的计算和观测水位变化趋势一致。由图 8可见，同时期的计算和实测平面流场基本一致，模型的拟合效果较好，符合国家标准^[25]，研究区地下水流数值模型的建立与实际相符。模型识别验证期地下水均衡为正均衡(表 4)，也符合研究区识别验证期末水位整体上升的实际情况。

4 地下水可开采量优化评价 4.1 优化管理模型

确定合理的地下水位对于研究区生态环境保护至关重要，而地下水位的高低直接受控于地下水的开采量及开采布局，为了实现地下水资源开发与生态环境协调发展，建立了基于水位限制的地下水优化管理模型和数值模拟模型进行耦合迭代计算，最终确定出羊庄岩溶水系统地下水的最大可开采量及最优开采布局。

将11个水源地的总开采量最大作为目标函数，把单个水源地开采量作为决策变量。考虑到各水源地均为生活和工业供水，需水量比较稳定，故将各水源地开采量均作为不随时间变化的常量进行优化计算，据此可建立优化管理模型。

目标函数为盆地内岩溶水开采量最大，即

式中：Q_总为11个水源地的总开采量；Q_i为第i个水源地的开采量; i为水源地编号。

约束条件如下。

1) 水位约束。各集中水源地在不同时刻的水位均不应低于该处最低约束水位，即

式中：H_i^t为第i个水源地在t时段水位标高；H_i⁰为第i个水源地的开采约束水位(表 5)。

2) 非负约束。决策变量应均为正值：

采用响应矩阵法将前面校正和验证后的地下水流模型与该优化管理模型(包括目标函数、约束条件)耦合起来，通过模型联合运算，便可确定出羊庄盆地地下水资源可开采量以及各个水源地优化开采量的空间分布。

4.2 优化管理模型条件确定

1) 优化模拟阶段降雨量

如何科学地预测未来的降水量，是地下水资源评价与合理开发的关键问题之一^[24]。此次评价根据羊庄气象站已有的降水观测数据，基于非平稳时间序列方法预测得到未来降水变化。可将月降水量的时间序列视为一维非平稳时间序列，其含有两种类型的成分：一种成分反映了序列随时间变化的平均规律，可用函数来逼近，该成分为决定性成分；另一成分是均值为0的平稳随机部分, 可用已知的线性模型来拟合。用公式^[26-27]可表示如下：

式中：T(t)为趋势项，反映降水量P(t)的变化趋势；Z(t)为周期项，反映P(t)的年际周期性变化, T(t)和Z(t)总体反映了时间序列变量变化中的确定性因素；剩余序列R(t)为随机项，反映了随机因素的影响，一般假设其为平稳随机序列。

笔者利用羊庄盆地1952—2015年的64年降水系列资料作为拟合数据，预测得到未来20年本区的降水可能值(图 9)，结果显示预测的20年平均降水量为707.74 mm/a，为一偏枯降水序列。

2) 开采约束条件

为实现羊庄盆地地下水环境持续好转的总目标，确定地下水资源优化配置的思路是，在目前分散生活用水和农业用水基本保持不变的前提下，以达到水源地总开采量小于盆地地下水总补给量，使得各水源地水位不低于约束水位，不诱发地质灾害为基本原则。

岩溶塌陷是近年来频发的主要地质灾害，其中地下水位降幅过大是岩溶地面塌陷的重要诱发因素之一。因此，为避免盆地内岩溶塌陷的发生，应该对各个水源地开采过程中地下水位最大降深或水位高程设定约束值。

2015年为羊庄盆地降水特枯年份，盆地地下水位历史最低，但整个盆地未发生岩溶塌陷，因此，以各水源地2015年的最低水位为基础：若其低于岩溶含水层顶板标高则可采用为约束水位；若其高于岩溶含水层顶板标高则采用岩溶含水层顶板标高为约束水位。最终根据各个水源地特枯年份(2015年)最低水位、岩溶含水层顶板深度以及含水层的开放性特征综合确定出水源地开采约束水位，见表 5。

4.3 优化管理模型求解及结果分析

使用Visual Modflow中的MODFLOW2000模块建立MODFLOW模型，经过识别与验证之后，从中提取GWM所需的MODFLOW数据文件，然后在GWM软件中输入优化计算所需要的决策变量、目标函数、约束条件以及优化方法，最后利用GWM所提供的序列线性规划法来求解该优化模型，从而得出羊庄盆地最大可开采量为5 730万m³/a，其中各水源地开采量优化配置方案见表 6。

根据优化后的开采方案，模拟预测了各水源地20年地下水位变化情况(图 10、11)，由模拟结果可知各个水源地岩溶水位均未出现持续下降的现象，且各个水源地开采水位最大深度亦均未超过拟定的埋深约束(表 5)，预测期20年研究区为正均衡(表 7)，说明本区岩溶水在含水层强大的调节作用下，5 730万m³/a可开采资源是有保障的。

5 结论与建议

1) 利用羊庄盆地的自然边界，对整个盆地进行建模，模型计算所得观测孔水头的变化趋势与实测水头基本一致。将地下水数值模拟模型与水资源优化模型进行耦合迭代，采用优化指导模拟、模拟验证优化的思路对开采方案进行合理调控，最终评价出盆地最大可开采量为5 730万m³/a，进而确定出各水源地的优化开采量。各个水源地岩溶水位均未出现持续下降的现象，且各个水源地开采水位最大深度亦均未超过拟定的埋深约束。

2) 鉴于大气降水的多少及其分配情况是影响地下水资源的最重要因素，此次基于非平稳时间序列的降水量预测是建立在对历史数据地分析之上，根据已有的降水数据预测而得到的未来降水变化，其摆脱了人为构造降雨序列的随意性，因此以该降水序列为依据模拟计算得到可开采量更为合理。

3) 由于资料有限，而且羊庄盆地尺度较大，此次羊庄盆地地下水数值模型的构建主要是基于等效多孔介质理论，该模型的局限性是未能精确地刻画局部含水层的高度非均质性。在未来的研究中，可以进一步调查该区岩溶裂隙发育程度以及发育主方向，并根据资料收集情况，综合选用双重介质模型或等效裂隙管道模型对局部重点地区进行小尺度模拟，以进一步丰富对该区域水文地质条件的认识。