地下水系统最优开采数值模型.pdf

第1 2 卷 第1 期 1 9 9 0年 2月 西安地质学院学报 J OURNAL OF XI AN COLLEGE OF GEOLOGY Vo 1 ． 1 2 No ． 1 Fe b． 1 9 9 0 地下水系统最优开采数值模型 仵 彦 卿 西安地质学院水工系 摄■本文以河南汝河 中段河各浅层地下永 系统为例，重点介绍 了一种地下水系 统 最 优 开采数值模型。该模型是一种地下水 减系统水均衡模型和线性规划模型的耦台 水 均 衡 模型采用强隐式有限差分法 S I P 求解，以解决技术函数及地下水 流系统势能 分 布 问 题’线性规划模型 以总开采置最大作为目标函数，以流量及技术函数组成 的响应矩 阵 作 为约束条件，甩单纯形法求解。该模型能比较快地解决多节点 薄层 潜水含水层中 地 下 水最优开采的规划问题t能提供 最优开采方案及最优开采量等信息 关■词水均衡模型J线性规划模型} 目标函数’技术函数’最优开采数值模型 引 言 地下水资源管理包括两方面涵义t一是水文地质工作者根据野外调查和试验所得资料， 运用系统分析方法，提供研究区地下水资源最恍开采方案，二是水文地质管理人员运用政治 的、法律的 、社会 的手 段实施最 优方案 。二者的 结台 才可实现对 地下水 资源的 有效 管理 。作 者提出的地下水系统最优开采数值模型解决的问题，应属于管理 问题的第一方面。运用该模 型研 究 了汝河 中段 河谷 浅层地 下水 系统最优开采 问题 ，通过 运营模 型提 出了研 究区未 来地下 水资源最优开发方案，即最优开采量、最优开采部位及开采方式。 浅层地下水系统概 述 研究区位于河南省郏县境内，属平顶山市管辖。由于平顶山市工业发展、矿 产 资 源 开 发，需水量增加，所 以开采该区地下水对城市建设尤为重要。 研究区面积1 8 1 x 1 0 m ，属半 湿润气候环境，大气降水多年平均7 0 0 ra m。地下水系统输 入 以大 气降 水为主 ，输 出以排 向汝河为主 。汝河纵 穿全 区，常年 有水 ，最 大 洪水期 流量 3 0 4 0 m ／ s ，最小 枯水期 流量5 ． 5 1 9 mS ／ s 。在未来集中开采条件下，汝河是主要的补给源。 含水层系统的物质分布主要是Q Q 。 河流冲积相砂 、砾卵石层，物质颗粒 粗 大，分 选 磨圆较好，厚度由西、北向南、东方 向逐渐变厚，最厚2 5 m，最薄5 m左右，平均约1 5 m。北 部为Q 冰碳物组成的垄岗，其底都为Q 钙质牯土或含卵砾石牯土，具有一定胶结， 形成稳定 1 9 8 9 一 O 3 0 6收 稿 维普资讯 1期 仟彦卿 地下东系统最优开采数值模型 隔水底板 ，南，西 ，东部 为人为边 界 ，其 物质 裹 1 研究区地下水系统参敛分布 组成 与外系统连续 。Ta b l e 1 Ds l t r i b u t i 。 “。 f p a r a me t 。 。o f 地 下 水 流 系统 属 开 放、 交替 积 极 的 水 动 力 兰 三 二 竺 竺 竺 系统 。上部 为潜水 自由面边界 ，北部与底部为 舒区号 I Ⅱ Ⅲ Ⅳ 零 流量 边界，西 北方 向为定水 头边界 开采 未 K 3 3 ． 2 ． 0 4 2 3 0 影响区 ，东部为定流量边界， 南部汝河作为定 。 - 。 _ 。 。 水 头边界。地 下水水 动力参数 k、 如 表1一 所 示 。 出 i 水 由于地下水 流积 极，形成 低矿 化度淡水， 藉 滤作用形成 的 HCOs - --Ca “ M 型水 ，是 比较 好的水源 。 综上 所述，研究 区地 下水 系统的特点姓 含水 层系统分布面广 ，物质 颗粒粗 ，孔隙介质通 畅 ，但 厚度较 薄。地 下水流系统可用一 确定性 、非稳 定流分布参 数型系统描述 见 图1。 b 图 1 研 究 区浅 层 地下 承 系统 圈 Fi g．Gr oundwat e r s y s t e m ma p of s hal 1 0 w i n s t u di e d a r e a ． 平面 目 b． 剖面 图 l 一潜水辱水位线I 2 一基岩探露I 3 一地质孔 ；4 一长观孔 I 5 一水文孔 t B 一隔水边界 ； 一定水头边界 ，8 一参数分区 透镜状粉钿砂 I l 0 一砂．砾 ，卵石含水层 I 1 1 m表土层 I 1 2 一粘土隔水罄 ，I 3 一流出边界1 4 一砾质粘土县 地下水流系统数学模型 0 ． 1 璩模型 建立 由上述地 下水 系 统特征 描述及 图1 示，研究 区地 下水 流系统数学模 型可 陈述为 维普资讯 舀安 地 学 琏学 扳 1 2 期 蛊 [ t 一 避】 【 h - q 嚣 卜N x ， ，t 。at ， y ∈Q 1 t ≥ 0 o c ≥ o hl F Ih 0 t ≥ O hI t ． 。 h 。 x， Y ∈Q k h 誊 一 c ≥ 0 或 ， 州 h～ B一Q 2 3 4 5 式中，在某共同的基准面上，h x ， y ， t 和 x ， Y ，h o x ， Y 分别为潜水面，隅水底板及初始 水 位标高 L ，r 、rz 、 和 Q分 别为定水头边界 、零流量边界 、定流量边界 和渗流区 域， k 和 分别为 渗透 系 数 LT t 和给水 度 无量纲 j Qj 为j t 单元开采量 L 。 TI 1 j 1 ， 2 ⋯M ， M为开采单元数 I N j x， Y 。 t 为垂 向天然补给强度 L T- 1 ，q x ， y 为单宽 流 出量 I T t }8为D r a c 的＆ 函数 L }n 为 内法 线方程 ； B为过 水断面宽度 ， 方 程 1 为非线性微分方程式，在求解 中比较麻烦。尽管如此 我们仍然可以采用 数值逼近方法使其 分时段线性化， 在每个时段内取平均厚度作为该时段内厚度值， 即T k [ h x ， Y 。 t 一 x ， y ] ≈K 丽 [ M 2 x ， y， A t I 十M3 x， y， At a ] ／ 2 ，M2 x ， ly ， △ t n 2 1 h x ， Y ， △ t n 一 Ⅵ x ， y ，M。 x ， y ， △ t h x ， y t a 一 x ， Y 。这样处理，其近似程度高，并保证计 算 中TK姐 常数 ，于是式 1可转变 为 3 ． 2模型 分解 为 了计算 技术 函数的需要，必须 把式 2～ 6所 组成 的定解 问题 分解 成如 下两个 定 解 问题 Ⅱ 基 T 3 、 T 3dh I ， x , y , t 一 量 Q X - X ,y ￡ ≥ 0 x ， Y ∈Q h l l t ． 0 h 。 X ， y ∈Q h l 1 r I h 0 t ≥ O 。 t ≥ o y X X ＆ Q M∑ 一 ， n 一 l、 a Q ∈ ， 嘧 Q M∑ H y Z N } 二、 a ，●＼ d } 、●， 扭 ， ● ＼ 弧 维普资讯 1期 仵彦卿 地下东系统最优开采数值模型 式中， Q 6 x x r ， y - y ， 是把r s 边界上单宽流量q x ， Y 转换成零流量 F 1 边界而 其 有 垂 向同样量 的单元 抽水量 。 f 鑫 T T 努 ～ Q j1 s cx 一 ， j1 等 I{ h Ⅱ I ．。 。 x ， y ∈ Q ≥ 。 ’ ∈ Q I h Ⅱ 】 r - 0 t ≥ 0 l o t ≥ o 模型 Ⅲ 是一个齐次初始条件。齐次边界条件下， 纯粹 由人工抽水引起的流场分布， 可 以证明h l 为负值，且是降深值。所以，在计算导水系数T时，应特剐注意含水层厚度，此 时 Tk [ h a x ， Y ～ x ， Y h Ⅱ x ， y ， t ] ，其它方法同上。 ●． 3 数 值 解 式 6可用 艟式 有 限差分格式 ； B i ， j h ’ D i ， h t El ， ．h ‘ ”F i ， - h l Pl ， J h ．⋯ - 1 A ’ 0 来 逼近 。 式 中， _ ， l Ki ， i l ， 丽 ’ J _ 1 ， 2 T； ， J 一 1 ， 2； Dj ， j K i I ， ， i 蚵 i I ， i ， I Tl 一 1 ， 2 ， {， Fi ， j Tj 1 ， ， J， Pi ， j T‘ ， j 1 I 2 ， EI ， 一B i , j -Di , i -Ft j -P{ ， ， ， ， ， Tl ， j - 1 ， l 2 Ti ， j T⋯ 一 I／ T ， j T。 ， _ - 1 ， Ti 一1 ， 2 ， i 2 Ti ， J Ti ． I ， I ／ T ， I Ti 一 1 ， J ， Ti ， j l ， 2 T； ， j Ti ， 』 l ／ Ti ， J T￡ ， J I ， T| 1 ， j 2 T⋯ T r ． 1 ， I ／ T⋯ T⋯ ， j ， 7 厶 x 为 空 间 步 长 L l A Q j’ j _’ j 警h t l △ t n 为 n 时 刻 时 间 步 长 ， Q ⋯ 为 i ， j 点垂向量 T 。 t ， i 表示节点号， 12 为时段计数器。 式 7写成 矩阵形式 为； E B3 { h} { A } 0 8 或I 维普资讯 68 El FJ Di Ei 0 D。 B‘ 0 西安 地面 学 藐学 报 0 Pl F2 0 P E3 Fs 0 P8 D‘E．F． 0 P‘ hl h 2 h 3 ● l 2 卷 Al A2 As ● A N 上 式系数矩 阵[ B J 为一五对 角型矩 阵，用Lu分解法计算工 作量 甚大 对 多节点 问 题 而 盲 ，若给系数矩 阵中增加一小矩阵E C ] 为两对角型矩阵 ，形成下式t [ B]EC] CN D 采用 下列式造代 计算 ； [ B ] E C3 { h ’ } E B ] E C] { h } 一 E B q { h K J } { A} 9 式 中， k为 迭代计数器 。 该算 法称 为强 隐式迭代法 SI P，具体方法详见 薛 禹群编 水 文地质 学数 值法 。 最优化开采模型 水 源地勘察的最 终d的是 在技术经济条 件下提 供最大允 许开采量 撂此 最优开 采模 型的 决策 变量 为单元开采 量Q_ j t 1 ， 2 ⋯M ，M为开采 单元总 数。 目标 函数是开采 量 最 ， 即 l 一 一 一 维普资讯 1期 { 午彦卿 求系 统 最扰 } } 采 数 值棋 型 B g l M a x Q M a x ∑q ． 1 0 l h I 1 I 约束条件 Ⅳ M N 1 1 3 M l,h , A t Q j 。 ≤啦 一h 1 1 I． M 、 I 1 蓍 - ≤ D I 1 2 O ≤Qj ≤Q。 j 1 ， 2 ， ⋯M 式中，QIl为j 。 单元抽水量 L 。 T }D为总开采量 上 限 T 一 ，取 2 ． 3 1 5 m3 ／ s ； 为 单元开采 量上 限 I J 日 T一，取 0 ． 1 1 5 7 m ／ s ，h 为模型 Ⅱ求 得预报期Mt 观 测点水头 值 LI h 表示执行计划作用量后，在观测点M- 处 【 起 的 允 许 水 头 值 L，该 值 的 确 定要考虑下列几个因素；①地下水位不能低于含水层厚度1 ／ 3 ，②采用潜水泵 抽 水，地 下水 位 不能低于泵 的扬 程}◎要考虑 开采时总费用 I ④ 要考虑补给条 件，补 给 良 好 一地 段，水I 头 h 要求大点 ，否则小点 。 ， p M ， j ， At n 。 为技术 函数。其定 义为对 所研 究的地下水水流 系统， 在齐 次初始 及边 界 h x ， y ， t l o ， h x ， y ， t l r - o ， ‘ 盖 i 0 条 件下 ， △ t n 。 时 刻 ， j - 单 元 施 加 单位抽水量 这里取1 0 0 0 0 r n 。 ／ d ，M。 单元的降深即为B M。 ， j ， △ t n ，即称为技术函数。 求 解 术函数的方法用数值法，采用模型 I求解，其步骤为t 第一步先 调试 初始 场，使其 达封水量 均衡 这里h x， Y， t j 一 0。 第二 步用SI P法求 解 I ，在第 一时段 At 一 个特定井上拖加一个 单位抽水量 ，在 第一时段之后不施 加抽水，求 得△t 时刻 技术 函数为 p 1 ， 1 ， At 。 日 2 ， 1 。 At B M。 ， 1 ， At 1 8 1 ． 2 ， At ， p 2 ， 2 ， At ⋯ p 1 ， j ， At 1 ⋯ 8 2 ， j 。 ， At 日 h 2 ， A t 。 ⋯ p M j 。 ， A t 第三步方 法同上 ，在第二时段 △ t z 一个持定井上分别 施加一个单 位抽水 量 ，在 其 它时段不施加 抽水 ，求得 △ t z 时刻技术函 数为， 依 次类推 可确定口 M- ， j ， At w 为 日 Ml ， j l ， △t I 日 1 ， 2 ， A t ⋯13 1 ， j ， A t 、 p 2 ， 2 ， A t z ⋯13 2 ， j ， A t 2 l p M 2 ， A t ⋯fl M j 。 ， A t 』 日 1 ， 1 ， A t N ． 日 2 ， 1 ， At N ， fl Mh1 ， At j p 1 2 ， A t N 1 ⋯ 13 1 ， j ， A t 1 13 2 ， ， △ 邮 △ t N 1 I p M i 2 ， △ t 。 ⋯ p p ， j 。 ， △ t ． j ＼ ●● ●● ● ● ， ， _ ● _- ●●● ●●● ●●● -●1 k △ M 日日日 △ M 维普资讯 7 O 西 安地 质 学 院学 撮 第 四步 t在整个开采 期 N- At 内，技术函数为 ； N1 ∑ 口 】 ， 1 。 A t n n t I N1 ∑ fl 2 ， 1 ， A t n n 2 l N1 ∑ B M ， 1 ， A t n 。 Ⅱ2 -1 N1 N B 3 ， 2 ， At n 。 ⋯ 日 2 ， j - ， At a 。 nt I I 1 5 1 N1 Nl B Ml ， 2 ， A t n 。 ⋯ 日 Ml ， j - ， △ t n 。 1 3 把 式 1 3代入 式 儿 后 ，最优开 采模型 即可建 立 ，该 摸型在数 学规 划 中称 为线性规 趔模型，运用数学规划中解线性规划的单纯形法，编制程序计算。整个地下水系统最优开果 数值模型犏成一个主程序与 四个子程序 用F OR TR AN 言编程 ，在UV 6 8 微机上实现， 其计算流程如图2 t ● 输人剖面单元 几何 量 及 M K． ． h h ． 及 边 界条件 调用子程序 l 调用 予 程 序 2 输 入 目标 函数 及约 束 条件 兜 午 肄 深值 等 调用子程序3 调用子 程 序 4 调用hi ． ML ． j 】 N 】 At 及 Qj 值 子程序 1 用 S l P法计算1 1 得hⅡ 子 程序 2 用 S I P法计算 fⅢ 得 口 M| ． j l N】 At 子 程 序 3 用 单纯形法计算总开采量 Q及 qj ． j 1 ． 2 。 ⋯ ⋯ M。 打 印Qi 及 Q 子程序4 用 S I P法 计 算 c 1 得h 打印预报 水位hi j L END 图 2 计算 机程 序 框图 Fi g． 2 Fr a me l r gur e 。 f compt t t er pr og r a m n 1 瞰 吼 ∑州 △ 2 ∑ 【 j △ M 日 ∑ 维普资讯 1期 许彦 卿 t地 下 盎 系 统最 优开 采数 值模 型 计算成果分析 计算时采用 正方 网格 剖分，步长 △ X1 0 0 0 m，割分 单元 1 8 1 个 ，节点2 2 6 - ，时 间 步 长 选用 两种，一是 等步 长 At 9 1 ． 2 5 d，是 变步 长 At n 1 ． 2 5 1 ． 5At n _ 1，预报 期 1 4 年 。初始流坜选用 枯水 季节 1 9 8 6 年1 2 月 1 日医内观 测孔 实测 数据 ，采 用趋 势面加残 差 的 方祛模 拟初始 场，其结果见 图3 。 冈 3 初抬条停下拳头分布 1 9 8 6 年1 2 月 1 日 Fi g． 3 He ad dl s t r i bu t i on und g r i ni t i al c ondi t i O R1 9 8 6． 1 2． 1 水 文地质参 数的选取 ，运用区 内长观 孔1 4 个月观 澜 l 序 列，用常 规方法 EE h 计一h 实 为最小 反演 确定参 数 k ， 其 五个区的参 数值见 表1 及 图】 。 从优 化结果看，1 9 个优 化区均位于河边 且含水 层厚度 比较大 的地方 ，有利 于最大 限度夺 取激发的河水补 给量，与 实际 人们 的认识吻台 ，如 图4 所示 。 计算结果表明，该区最大允许开采量为1 ． 8 5 m ／ s ，其各优化单元开采量见表2 。 在最优开采条件下，预报了未来1 4 年枯季 2 0 0 0 年1 2 B 1 日地下水水头分布 见图4 ． 显见 ，地 下水水 位在西 部变 化不大 ，仅在开 采区的东 部形成 降落漏斗， 获得较大 的激 发河水 补给量。其值为1 ． 0 1 5 m ／ s ，占总补给量5 4 ． 9 ，河流补给保证程度达l O 0 。从流入流出及 变化量看，区内水量是均衡的 Q总 补2 ． 1 8 m ／ s s Q总 排2 ． 1 8 1 m 。 ／ s 。从降落漏斗中心选 取 1 、2 号观 测孔 见 图4，其 水 头随 时问变化 瞄线如图5 、图6 所示 。可见 ，地下水 位 最 大 降深4 ． 6 7 m，1 号节点最大降深 为4 ． 1 4 8 m，6 年基本 稳定 。 维普资讯 西安 翘质 学 院学 报 1 2 ； g 图 4 预报优化开采条R - T的水头分布 2 0 O 睥 l 2 月1 日 Fi g． 4 H e a d di s t r i but i on unde r o pt i mal de v el opme nt2 0 0 0． 1 2． 1 囊 2 优化簟元开采■曩辟霖 Ta bl e 2 Opt i mal yi e l d a nd dr awdown i n op t i mal e 1 e me nt h m} 9 2 ． 8 2 3 8 6 ． 8 2 3 7 3 0 l | l 6 O 2 1 9 0 2 ∞蛹O ‘ 3 H 5 t I o 图 5 1 号 观测 孔h ～ t 关 系 曲 线 Fi g． 5 h ～ t c t l r ve of No． 1 no de t B 维普资讯 1期 仟彦卿地 下水系统最优开果数值模型 7 3 m l ∞ ． 捌l 7 95 ． ∞ 7 O 1 3 0 1 0 0 图 F i g． 6 21 90 2 92 O 3 65 0 43 8 0 51 t O 6 2 号 观涮 孔h ～ t 关 系曲 线 h～ t c ur e o f N o． 2 n ode 结 论 t I B 1 ． 地下水系统最优开采数值模型的运用，可以达到经济合理地开采地下水资壤，避免由 于盲 目开采造成 的一 系列环境危害。 2 ． 对 于多节点 、非线性 比较 明显的 情况 如 大面积 ，薄层潜水 问题 ，采用 文 中计算 方 扶速 度快，且选择方案 灵活 方格 网剖分 选用 时间步长最好用变步 长。 3 ． 应 用该 模型计算，得 出研究 区最 优开采量1 ． 8 5 m ／ s 。 感谢 河南第三水文 队汝河组的大力 掷助和李 俊亭 副教授 在计 算方 面的具体 指导 参考文献 I 雅 贝尔著 许稻铭等译 ．地下亦亦力学 ．j 京-地质出麓社 9 8 THE NUM ERI CAL M 0DEL 0F 0PTI M AL DEVEL0PMENT 0F THE GR0UNDW ATER S YSTEM W u Yonqi * o Xi a Co l l e g e of Ge ol o g y Ab ； t r ac t Gi v i g an e xa mpl e a b ou t t he g r o un dwat e r s y s t e m o f t h e s h a l l ow of t h e H e na n Ruh e mi d dl e r e ac h． t hi s pa pe r h as ma i nl y i nt r o d u c e d a n u m e r l c a l mo d e 1 of op t i ma l d e e 1 op me n t of t h e g r o u n dwa t e r s ys t e m ． Th i s mod e l i s mad e u p of t h e wat e r b a1 a nt e mod e 1 of t h e g r o un d wa t e r f l o w s y s t e m a nd t h e l i ne a r pr og r a mmi ng mo d e1 ． Th e wa t e r b al a nc e m o d e l i s s o l v e d wi t h t h e SI P Fi n i t e Di f f e r e nc e me t h od s o a s t o co mpu t e t h e t e c hn ol o g i c a l f u nc t i o ne an d t h e po t e nt i a l e ne r g y d l t r i b u t i o n o f t h e g r ou n dwa t e r s ys t e m Th e l i n e a r 维普资讯 74 p r og r a mmi ng mode lo pt i ma l M d e v e l o pme nt mod e 1 t ak e s t h e t o t a 1 yi e I d Ma x ∑ Qj 1 a s a n o b j e e t f u n c t i o n a n d t h e R e s p o n s e Ma t r i x ma d e o f t h e j 1 - t e c hn ol og i c al f u nc t i o n a nd t h e f l ow a mo u nt a s a c o ns t r a i n t ． Th i s mod e 1 i s s o l ve d by t h e Si mp l e x me t h o d ． 一 W i t h t h i s mo de l t h e pr o g r a mmi ng pr o bl e m a bo u t t h e o pt i ma l yi e l d of t he g r o u nd wa t e r f r om t h e mul t i no d al po i n t s a nd t he t h i n p hr e a t i c a q u i f e r c a n he s ol v e d s o o ne r ， an d t h e i nf o r ma t i o n a bo u t t he op t i ma l yi e l d pr o g r a mm a n d t h e o pt i ma l yi e l d a mo u nt wi l l b e p r oY i d e d ． K a y wo r d s wa t e r b a l a n e e mo d e l ，l i n e a r p r o g r a mmi n g mo d e l ，o b j e c t f u n c t i o n， co ns t r a i nt ， t e c h no l o g i c al f u ne t i O n， o pt i r fi a l d ev e l op me nt nu me r i c al mo de l ， g r o u n dwat e r s ys t e m 维普资讯