迭代注水算法

‘壹’ 二维承压水非稳定流模型遗传反演方法

11.2.1 数学模型^［49］

含水层参数识别方法

其中^[29]：H 为地下水水头[L]；S 为贮水系数[量纲一]；T=KM 为导水系数[L²/T]；K 为含水层渗透系数[L/T]，二秩对称张量；M 为含水层厚度[L]；（x，y）为笛卡儿坐标系的坐标；t 为时间；G 为不含边界的研究区域；Γ1 为第一类边界；Γ2 为第二类边界；Γ=Γ1+Γ2 为边界；-G=G+Γ为包括边界的研究区域；n 为边界单位外法向量；H₀ （x，y）为初始水头分布；H₁ （x，y，t）为第一类边界上的水头分布；q（x，y，t）为第二类边界单位宽度上的流量分布；E 为垂向补给强度。

含水层参数识别方法

其中，N_W抽水井（或注水井）数目；Q_j为j号井中的出水量。

11.2.2 模型的解

（1）化边值问题（11-1）为变分问题

若固定 t，并使 H（x，y，t）∈C² 则对任意η（x，y）∈则有：

含水层参数识别方法

利用Green公式，并注意边界条件得：

含水层参数识别方法

（2）剖分求解区域

在区域G的边界Γ=Γ₁+Γ₂上，取有限个点依次联成一闭多边形Γ_D。以此近似代替Γ，并以Γ_D围成的多边形区域G_D近似代替G，然后把G_D剖分为N_e个三角形之和，三角形的顶点为节点，节点编号i=1，2，3，…，N_P；其中内部节点N₀个，第一边界节点N₁个，第二边界节点N₂个。第i号节点的坐标记为（x_i，y_i），三角形单元用Δβ（β=1，2，…，N_e）表示。

（3）构造基函数和线性形状函数

对以顶点为i，j，k的三角形Δβ上的任意点P（x，y）的面积坐标定义为：

含水层参数识别方法

式中，Δβ为三角形面积，Δi，Δj，Δk分别表示三角形pjk，pki，pij的面积（见图11-1）。

含水层参数识别方法

图11-1 三角形面积坐标示意图

含水层参数识别方法

节点i的基函数定义为：

含水层参数识别方法

其中，{Δβi}为以 i 节点为公共节点的三角单元的集合。

取形状函数为线性元，则在节点i上水头的分片线性函数可表示为

含水层参数识别方法

取在边界Γ₁上为0的任意分片光滑函数η为

含水层参数识别方法

式中，m=N₀+N₂，r_i为任意常数。

（4）建立有限元方程

将式（11-8）和式（11-9）代入（11-4）式并考虑r_i的任意性，得有限元方程：

含水层参数识别方法

考虑初始条件，用全隐式求解常微分方程组（11-3），可得：

含水层参数识别方法

式中，=H（x_i，y_i，t_n），Δt=t_n+1-t_n，m=N₀+N₂，m 为未知节点的个数。其中：

含水层参数识别方法

其中，N_t表示以i和j为公共节点的单元个数，T_xxβ，T_yyβ，S_β为T_xxT_yy，S在Δβ内的值，l_ij为i，j两节点之间距离，q_ij为线段ij上的流量。

含水层参数识别方法

其中，N_t表示以j为公共节点的单元个数，Q_j表示j节点处的抽水流量，Q_β表示j节点相邻单元中（x_β，y_β）处的抽水流量，Q_j=0表示j节点处无抽水，Q_β=0表示单元β中无抽水。

方程组最终可以写为如下矩阵方程：

含水层参数识别方法

A为大型稀疏对称矩阵，方程组（11-12）可用一维压缩存储技术直接求解，也可用高斯-赛得尔迭代方法求解，按不同的时间步长解出每一预定时刻的每一节点的水位。

11.2.3 反演方法简介

二维承压水非稳定流模型的反演方法，在第1章已经进行详细的论述，许多方法也比较成熟。在实际的水文地质计算中，最常用的是试估-校正法，但是该方法无收敛判别准则，很难达到最优识别，工作量比较大。使用该方法，结果的可靠性和花费时间的多少取决于调参者的经验和技巧。

11.2.4 遗传反演方法

遗传反演方法是一种最优化方法，是将水文地质参数识别归结为求极值问题。即求水文地质参数使得误差评价函数达到最小。

设共有n个水文地质参数，用符号p₁，p₂，…，p_n来代表。同时假设在j号观测点上i 时刻的计算水头为（t_i），实测水头为（t_i ），比较的观测点总数为 N 个，比较的时间段为 M 个。通常用平均误差绝对值及平均误差平方和来表示拟合的程度。称为评价函数（也称目标函数）E。显然，E 是所给出的参数值的函数。评价函数的表达式如下：

含水层参数识别方法

遗传反演方法是将反演求参问题转化为优化问题，以简单的遗传算法为基础，并对简单的遗传算法进行改进，有效地求解优化问题，从而反演求参的方法。

在前面我们已使用遗传反演方法对一维地下水系统的反演问题进行了求解，本节讨论的问题和前面讨论的问题没有本质的区别。只不过问题更复杂，所求参数由于参数分区的存在，不再是一组参数，而是多组参数。其反演过程前面已有所述，我们仍然对8种遗传算法进行讨论和比较。

a.简单遗传算法（SGA）；

b.优体克隆遗传算法（The Best Chromosome Clone GA=BCC-GA）；

c.优体克隆+子体优生（Younger Generation Chromosomes Prepotency）的遗传算法（BCC-YGCP-GA）；

d.优体克隆+子体优生+多代调环（Multi-Generations Adjusting Environment=MAE）的遗传算法（BCC-YGCP-MGAE-GA）；

e.多代调环的遗传算法（MGAE-GA）；

f.优体克隆+多代调环的遗传算法（BCC-MGAE-GA）；

g.子体优生的遗传算法（YGCP-GA）；

h.子体优生+多代调环的遗传算法（YGCP-MGAE-GA）。

例1.假设承压含水层区域是一边长为a的正方形，东西边界为定水头边界，水头为H₁，南北边界为隔水边界，区域中心有一抽水井以流量Q抽水，承压含水层的导水系数为T。非稳定流定解问题如下：

含水层参数识别方法

此定解问题的解析解由Chan^［39］，Mullineux和Reed给出。其中：（x₀，y₀）为抽水井的坐标。在计算时，正方形的边长a为1200 m，边界AB和DC为隔水边界，AD和BC为定水头边界，H₁=100 m。计算剖分图见图11-2，剖分三角形单元数为312个，节点数为181个，其中内部节点数为133个，一类边界节点数为26个，二类边界节点为22个。计算时分为三个参数分区，见图11-3（表11-1），在区域中心P₁点设置一抽水孔，抽水流量Q=1000 m³/d，T的单位为m²/d，S为量纲一变量。各参数分区的参数见表7-1，H的单位为m，含水层的顶底板高程分别为0和50 m，含水层厚度为50 m。在三个参数分区分别设置一个观测孔OBS1、OBS2、OBS3，其位置见图11-3，其水位观测值见表11-2。根据水位观测值反演水文地质参数。

图11-2 地下水水流模型有限元剖分图

图11-3 模型参数分区图

表11-1 模型各分区参数表

表11-2 观测孔水位观测资料表

在用遗传算法反演参数时，在所有方案中均采用：导数系数 T₁，T₂，T₃ 的初始取值区间为（0，1000.0 m²/h），贮水系数 S₁，S₂，S₃ 的初始取值区间为（0，0.1），遗传代数Num Gen=1000，种群数 Pop Size=50，交叉概率 P_c=0.7，变异概率 P_m=0.3，评价函数中的alph=0.05。所有与多代调环有关的计算方案，其代数选择均为10。前500代每隔10代对 T 和S 设置一次取值区间，取值区间设为（0.75V_best，1.25V_best），V 为变量可代表T 或S，V _best为这10代中最优的染色体，500代后开始进行区间压缩技术，其压缩方法见前章所述，压缩系数为0.4。本问题为理想模型问题，我们用8种遗传算法或其不同的组合进行计算其结果见图11-4和表11-3。

图11-4 遗传反演方法进化代数与目标函数计算结果图

表11-3 不同的遗传反演方法目标误差函数计算结果比较表

从计算结果可以看出与子体优生有关的遗传算法如BCC-YGCP-GA、BCC-YGCP-MGAE-GA、YGCP-GA和YGCP-MGAE-GA均达到了收敛，其反演结果均较好。其中最好的结果是YGCP-MGAE-GA，当达到700代时目标误差函数达到全局最小为0.000041，但这种方法不稳定很容易出现返祖现象。这种方法比别的方法表现好，可能包含随机因素在内。简单的遗传算法和优体克隆遗传算法在一千代时均不收敛，其反演计算结果与真值相差太远。在所有这些遗传算法中优体克隆+子体优生遗传算法收敛速度最快，当到100代时其目标误差函数为0.00129，其反演出的参数三个分区的渗透系数分别为：51.520，249.755，500.895；三个分区的贮水系数分别为0.000101，0.000207，0.000479。反演的结果已接近真值。在这里我们列出简单遗传算法、优体克隆+子体优生遗传算法、子体优生+多代调环遗传算法的计算结果以便比较。其结果见表11-4、表11-5、表11-6。

表11-4 简单遗传算法反演结果表

表11-5 优体克隆+子体优生遗传算法反演结果表

表11-6 子体优生+多代调环遗传算法反演结果表

‘贰’ 哪些职业看上去很“美”，前景也很好

职业的“美”与不美，每个人的定义是不一样的。就我而言，职业看上去很美，一定是指那些人人都羡慕的职业，这些职业不见得轻松，但一定有一些点让你非常向往，比如自由、高薪等等。所以我觉得看上去很美且前景比较好的职业有记者、演员、律师。

律师总是能西装革面，却能轻松赚到高薪。律师这个职业简直就是“社会精英”和“潇洒多金”的代名词，很多电视的男主都已经被律师职业承包了，所以说它“美”，一点儿也不过分。再者说，律师是越老越吃香，年龄越大，经历的案子多了，社会阅历也更丰富了，所以找你办案子的人就会越多。所以律师这个职业是很有前景的，谁不想动动嘴皮子就能拿一笔不少的报酬呢。

‘叁’ MIMO技术原理及应用的图书目录

第1章 绪论
1.1 MIMO系统的提出
1.2 MIMO系统的特征及研究进展
1.2.1 MIMO系统的主要特征
1.2.2 已取得的进展
1.3 存在的问题
参考文献
第2章 MIMO信道建模
2.1 无线信道建模的必要性
2.1.1 大尺度衰落及其典型模型
2.1.2 小尺度衰落及其典型模型
2.1.3 信道的一阶和二阶统计量
2.2 MIMO信道建模的研究现状
2.2.1 MIMO信道建模的必要性
2.2.2 从SISO信道到MIMO信道的演变
2.2.3 MIMO信道建模方法的分类
2.2.4 MIMO信道典型模型
2.3 MIMO信道建模两个实例研究
2.3.1 MIMO无线信道参数
2.3.2 MIMO信道的空间相关性
2.3.3 基于Kronecker的MIMO信道模型
2.3.4 单环及改进型单环MIMO信道模型
参考文献
第3章 衰落信道的容量
3.1 高斯信道下的信道容量
3.2 平坦衰落信道的容量
3.2.1 信道与系统模型
3.2.2 接收机知道信道状态信息，发射机知道信道分布
3.2.3 发射机与接收机均已知信道状态信息
3.2.4 分集接收机的容量
3.2.5 相关Nakagami信道分集接收机的容量
3.3 频率选择性衰落信道的容量
3.3.1 时不变频率选择信道
3.3.2 时变频率选择信道
参考文献
第4章 MIMO信道的容量
4.1 独立衰落下单用户MIMO系统的容量
4.2 信道系数固定时的MIMO系统容量
4.2.1 循环对称复高斯随机向量
4.2.2 通过互信息推导MIMO系统的容量
4.2.3 通过信道矩阵的奇异值推导MIMO系统的容量
4.3 信道系数随机变化时的MIMO系统容量
4.3.1 容量的定义
4.3.2 MIMO系统的各态历经容量
4.4 MIMO系统的容量实例及仿真分析
4.4.1 单输入单输出(SISO)系统的容量
4.4.2 SIMO系统的容量
4.4.3 MISO系统的容量
4.4.4 两种典型的MIMO系统容量
4.5 相关衰落下单用户M1MO系统容量
4.5.1 接收机能准确估计信道，发射机不能估计信道
4.5.2 接收机和发射机均不能估计信道
4.5.3 频率选择性衰落相关信道下MIMO=OFDM系统容量
4.6 多用户：MIMO系统容量分析
4.6.1 MIMOMAC系统
4.6.2 MIMOBC系统
4.6.3 MIMO-MAC和MIMOBC的对偶性
4.6.4 迭代注水算法
4.7 基于训练序列估计的MIMO系统容量
4.7.1 基于训练序列信道估计的MIMO系统模型
4.7.2 基于训练序列的信道估计值的推导证明
4.7.3 等效的系统模型
4.7.4 基于训练序列估计的信道容量
参考文献
第5章 分集技术
5.1 分集类型
5.2 分集增益与编码增益
5.3 接收分集系统模型
5.4.发射分集
5.4.1 发射机不知信道状态MISO
5.4.2 发射机已知信道状态：MISO
5.4.3 发射机已知信道状态：MIMO
5.5 矩分析方法在分集技术中的运用
参考文献
第6章 空时编码技术
6.1 空时编码技术基础
6.1.1 空时编码模型
6.1.2 空时编码的性能分析
6.2 空时编码设计准则
6.2.1 慢衰落瑞利信道的编码设计准则
6.2.2 快衰落瑞利信道的编码设计准则
6.3 空时编码的性能指标
6.4 空时编码的成对差错概率的准确估算
6.5 空时格形码性能分析
6.5.1 空时格形码的编码方案
6.5.2 空时格形码的译码方案
6.5.3 空时格形码的性能分析
6.6 基于正交设计的空时分组码
6.6.1 Alamouti发射分集方案
6.6.2 空时分组编码的正交设计
6.6.3 准正交空时编码的基本原理和设计准则
6.7 基于星座旋转的满分集的准正交空时编码
6.7.1 满分集的准正交空时编码设计
6.7.2 满分集的准正交空时编码的性能指标
6.8 空时编码器
6.8.1 空时信号的构建
6.8.2 空时码的性能
6.9 差分空时码
6.9.1 单天线系统中的差分空时码
6.9.2 MIMO系统中的差分空时码
参考文献
第7章 MIMO系统检测算法
7.1 单小区情况单用户MIMO系统模型
7.2 最大似然检测
7.3 线性检测算法
7.3.1 基于迫零准则
7.3.2 基于：MMSE准则
7.3.3 串行干扰抵消算法
7.4 非线性检测算法
7.4.1 QR分解算法
7.4.2 MMSE意义上的SQRD
7.5 结合格缩减技术的检测
7.5.1 基本原理
7.5.2 格缩减技术
7.5.3 格缩减辅助的检测算法
7.5.4 格缩减辅助的线性检测
7.5.5 格缩减辅助的BLAST非线性检测
7.6 球形译码算法(SDA)
7.6.1 FP算法
7.6.2 VB算法
7.6.3 SE-VB算法
7.6.4 自动球形译码算法
7.6.5 各种改进版本的k-bestSDA
7.7 Q1w算法
7.8 半定松弛算法
7.8.1 关于松弛的基本概念
7.8.2 半定松弛最大似然检测
7.9 分枝定界算法
7.10 堆栈算法
7.11 智能检测算法
7.11.1 禁忌搜索检测
7.11.2 粒子群优化
7.12 蒙特卡罗统计等算法
参考文献
第8章 MIMO中继信道
8.1 协同通信
8.1.1 协同MIMO技术
8.1.2 协同中继传输
8.1.3 用户协同传输
8.1.4 协同通信技术特征
8.2 加性高斯信道协同无线信道容量
8.2.1 三节点中继信道模型
8.2.2 半双工协同中继方法
8.2.3 半双工解码前向中继
8.2.4 半双工放大前向中继
8.2.5 半双工选择性中继
8.2.6 半双工增量中继
8.3 多节点高斯协同中继信道
8.4 衰落信道.MIMO协同中继系统容量
8.4.1 传统MIMO信道容量
8.4.2 MIMO协同中继系统容量
8.5 协同中继系统的功率分配
8.5.1 中继链路系统模型
8.5.2 中断概率相等功率分配策略
8.5.3 DF中继链路功率分配策略
8.5.4 AF中继链路功率分配策略
8.5.5 仿真分析
8.5.6 MIMO协同中继系统的功率分配
8.5.7 仿真分析
8.6 协同功率分配
8.6.1 三节点两跳中继网络
8.6.2 多节点两跳中继网络
参考文献
第9章 MIMO.OFDM系统
9.1 OFDM系统基本概念
9.2 OFDM的系统结构框图
9.2.1 OFDM主要功能模块
9.2.2 串并变换
9.2.3 子载波调制
9.2.4 DFT的实现
9.2.5 保护间隔和循环前缀
9.2.6 OFDM系统的缺点
9.3 基于IEEE802.16的WiMAX系统
9.3.1 IEEE802.16无线接入标准
9.3.2 WiMAX论坛
9.3.3 物理层关键技术
9.3.4 IEEE802.16物理层简单介绍
9.3.5 IEEE802.16e的网络结构
9.4 IEEE802.11无线局域网标准
9.5 LTE系统简介
参考文献
第10章 MIMO天线设计
10.1 概述
10.2 MIMO多天线与传统天线设计的比较
10.3 MIMO天线设计基础
10.3.1 MIMO天线单元设计要求
10.3.2 设计思想
10.4 天线设计准则
10.5 MIMO移动台天线设计
10.6 MIMO基站天线设计案例
10.7 多模式天线在MIMO系统中的应用
10.7.1 同轴波导馈电的双锥天线
10.7.2 自补偿阿基米德四臂
螺旋天线
参考文献

‘肆’ 注册岩土工程师考试科目

基础考试为闭卷考试：

上午段主要测试考生对基础科学的掌握程度，设120道单选题，每题1分，分11个科目：高等数学、普通物理、普通化学、理论力学、材料力学、流体力学、电工电子技术、信号与信息技术、计算机应用基础、工程经济、法律法规。

下午段主要测试考生对岩土工程直接有关专业理论知识的掌握程度，设60道题，每题2分，分8个科目：土木工程材料、工程测量、土木工程施工与管理、工程地质、结构力学、结构设计、土力学与基础工程、岩体力学与岩体工程。

专业考试的专业范围包括：岩土工程勘察、浅基础、深基础、地基处理、边坡和基坑、特殊土和不良地质、建筑工程抗震、地基检测。专业基础考试上午和下午各100分；

专业案例考试，上午段和下午段各有30个题目，从这30个题目中选择25个进行考试，上午下午各50分，共计100分。

(4)迭代注水算法扩展阅读：

报考资格

凡中华人民共和国公民，遵守国家法律、法规，恪守职业道德，并具备相应专业教育和职业实践条件者，均可申请参加注册土木工程师(岩土)执业资格考试。

(一)具备以下条件之一者，可申请参加基础考试：

1、取得本专业(指勘查技术与工程、土木工程、水利水电工程、港口航道与海岸工程专业，下同)或相近专业(指地质勘探、环境工程、工程力学专业，下同)大学本科及以上学历或学位。

2、取得本专业或相近专业大学专科学历，从事岩土工程专业工作满1年。

3、取得其他工科专业大学本科及以上学历或学位，从事岩土工程专业工作满1年。

(二)基础考试合格，并具备以下条件之一者，可申请参加专业考试：

1、取得本专业博士学位，累计从事岩土工程专业工作满2年;或取得相近专业博士学位，累计从事岩土工程专业工作满3年。

2、取得本专业硕士学位，累计从事岩土工程专业工作满3年;或取得相近专业硕士学位，累计从事岩土工程专业工作满4年。

3、取得本专业双学士学位或研究生班毕业，累计从事岩土工程专业工作满4年;或取得相近专业双学士学位或研究生班毕业，累计从事岩土工程专业工作满5年。

4、取得本专业大学本科学历，累计从事岩土工程专业工作满5年;或取得相近专业大学本科学历，累计从事岩土工程专业工作满6年。

5、取得本专业大学专科学历，累计从事岩土工程专业工作满6年;或取得相近专业大学专科学历，累计从事岩土工程专业工作满7年。

6、取得其他工科专业大学本科及以上学历或学位，累计从事岩土工程专业工作满8年。

(三)符合下列条件之一者，可免基础考试，只需参加专业考试：

1、1991年及以前，取得本专业硕士及以上学位，累计从事岩土工程专业工作满6年;或取得相近专业硕位士及以上学，累计从事岩土工程专业工作满7年。

2、1991年及以前，取得本专业双学士学位或研究生班毕业，累计从事岩土工程专业工作满7年;或取得相近专业双学士学位或研究生班毕业，累计从事岩土工程专业工作满8年。

3、1989年及以前，取得本专业大学本科学历，累计从事岩土工程专业工作满8年;或取得相近专业大学本科学历，累计从事岩土工程专业工作满9年。

4、1987年及以前，取得本专业大学专科学历，累计从事岩土工程专业工作满9年;或取得相近专业大学专科学历，累计从事岩土工程专业工作满10年。

5、1985年及以前，取得其他工科专业大学本科及以上学历或学位，累计从事岩土工程专业工作满12年。

6、1982年及以前，取得其他工科专业大学专科及以上学历，累计从事岩土工程专业工作满9年。

7、1977年及以前，取得本专业中专学历或1972年及以前取得相近专业中专学历，累计从事岩土工程专业工作满10年。

参考资料:网络----注册岩土工程师

‘伍’ 模型的特色及改进

一、混合井孔的模拟

为了简化成井工艺，增加出水量和降低成本，本区地下水开采多采用混合井（图5-2）。

混合井在地下水流模型中的刻画，目前大都采用由美国地调局1988年推出，并在国际上广泛应用的三维有限差分地下水流模型MODFLOW软件处理。

图5-2 混合抽水井

MODFLOW建议：多层井的流量必须以某种方式人为地分配给每一单层，……把井流量按每层的导水系数大小分配，即Q_i/Q_w=T_i/∑T（1988，1996，2000年）。其中Q_i和Q_w分别为第i层流量和总（井口）流量，T_i和∑T分别为第i层导水系数和总导水系数。

为了便于讨论，又不失其一般，我们以贯穿两个含水层的混合井孔（图5-2）为例进行讨论。如此，上式可表述为Q₁/Q₂=T₁/T₂。

MODFLOW对混合井的这种处理方法没有给出理论上任何的分析、说明，缺乏理论依据，应用中也与实际不符（陈崇希等，1998；Chen Chong-xi和Jiao J J，1999）。这是因为：

1）含水层导水系数对井孔流量的影响，不会如此简单。例如，混合井附近岩性（渗透系数）发生变化，甚至混合井打在岩性透镜体上（图5-3），怎样影响各层流量的分配?含水层厚度发生变化又如何改变流量的分配?MODFLOW无法解答这些经常遇到的实际问题。

图5-3 混合井管贯穿岩性透镜体

2）含水层的参数影响混合井流量的分配，导水系数只是其中一个因素，含水层的弹性给水度（储水系数）也应该起作用。

3）井管的流量分配不仅与含水层的水文地质参数分布有关，还与外边界条件、井径（有效井径）、水泵吸水管的位置及其他抽水井的干扰等因素有关，MODFLOW仅依导水系数对混合井的流量进行预分配，未能考虑其他因素的作用。

4）含水层的导水系数一般不随时间变化，即导水系数比T₁/T₂不会变化，但流量比Q₁/Q₂却不是一个常量。

5）按MODFLOW的方法，一个混合抽水井的流量Q₁/Q₂始终是个常量。然而在模拟过程中，特别是预测时，周围随时可加入或关闭混合的或非混合的抽水井，在这种井群干扰下，原混合井的流量比还会保持常量吗?显然是不可能的。

6）从另一角度分析，混合观测孔是混合抽水井的特殊情况（Q_w=0），对于两层混合的观测孔，其孔中水位（混合水位）必介于两含水层水位之间，即混合观测孔对于其中一含水层（例如1含水层）起抽水作用（Q₁＞0），对于另一含水层（2含水层）起注水作用（Q₂＜0）。如此，Q₁/Q₂＜0。而两含水层的导水系数的比值肯定是正值，即T₁/T₂＞0。如此，两个比值怎能相等?有人认为：MODFLOW并没有说上式可用于混合观测孔。对此，很容易论证：当混合抽水井的抽水量足够小，足以保持混合井中的水位介于上下两含水层的水头之间，上述论证一样成立。

7）按各层导水系数T的比例来预先人为划分各层的流量。这是缺乏理论依据的，因为这种做法要求各分层的有效井径相等和井壁处的水力坡度上下处处相等，这两个条件不可能人为控制，预先也不得而知。

显然，用导水系数的比值预先给定各分层的流量是不妥的。实质上，这种方法不是模拟，而是“处理”，一种未考虑机理的“处理”。

“防止地下水模拟失真，提高仿真性”是水文地质模拟工作者的核心任务，而地下水流系统中普遍存在的混合井，又是当前国内外模拟失真的主要问题之一，应当引起我们的足够重视。

研究区的混合井，既有抽水井，也有观测孔。三维流场中，即使属于均质含水层（无需多层含水系统），常规（理论上非点状滤管）抽水井和观测孔都属于混合井孔，因为滤管中不同深度处的水头是不相等的，因此滤管中的水要发生垂直流动，即使在不抽水的观测孔中也一样。这就是地下水流对混合井孔响应的本质所在。

本研究采用“渗流-管流耦合模型”（陈崇希等，1992，1996）来刻画混合井孔，以“渗流”刻画地下水的运动，“管流”刻画井孔中的水流，以解决混合井孔的模拟问题，大大地提高了模型的仿真性。“渗流-管流耦合模型”的基本思路如下（图5-4）：

图5-4 计算框图

1）将产生管流的混合井视为渗透系数很大的“圆柱形透镜体”，含水层通过这个具有很大导水性的圆柱体（井筒）强烈地交换水量，那么混合抽水问题就可以视为一个特殊的“越流系统”。

2）求越流系统“圆柱形透镜体”的“渗透系数”。需要强调的是，“圆柱形透镜体”的“渗透系数”不是常量，而是随流态的变化而变化。

具体的计算方法如下：

由流体力学知识可知，当管流呈层流状态时，其水头损失可依Darcy-Weisbach方程计算：

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：ΔH为水头损失；f为摩擦系数；l为管长，m；d为管内直径，m；u为管内平均流速，m/d；g为重力加速度，m/d²。

当管流为层流时

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：Re为雷诺数。

河西走廊疏勒河流域地下水资源合理开发利用调查评价

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：ν为流体的运动黏度；μ为流体的动力黏度；ρ为流体密度。

将式（5-2）至（5-4）代入式（5-1），得

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：J为水力坡度；γ为流体重度。

将上式写成渗透流速的形式，对于管流，其空隙率n=1，则渗透流速v=nu=u，故有

河西走廊疏勒河流域地下水资源合理开发利用调查评价

将此式与达西定律v=KJ对比，可得出层流状态下管流的等效渗透系数K_L的表达式

河西走廊疏勒河流域地下水资源合理开发利用调查评价

此式由陈崇希（1966）和J.Bear（1972）分别获得。

根据流体力学知识可知，当雷诺数大于100000时，摩擦系数f与雷诺数Re无关，而取决于管内壁的相对粗糙度（），这时的水头损失ΔH与流速的平方u²成正比。在3000＜Re＜100000的范围内，其中存在f=，即ΔH∝u^1.7⁵的区段。在紊流条件下，除了上述两个区段外，还存在两个过渡区，共分为4个区段。可见，紊流是个比较复杂的问题。为了解决此问题，这里提出等效渗透系数K_N的概念。

当管流呈紊流状态时，式（5-1）可改写为

河西走廊疏勒河流域地下水资源合理开发利用调查评价

如果定义紊流状态下管流的等效渗透系数

河西走廊疏勒河流域地下水资源合理开发利用调查评价

则

河西走廊疏勒河流域地下水资源合理开发利用调查评价

该式也具有达西定律的形式。

本来，紊流的运动规律有别于达西定律，而且不同流态区具有不同的形式。引入等效渗透系数后，将5个流态分区（1个层流区和4个紊流区）的运动规律统一为达西定律形式，且与地下水渗流定律一致。即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

其中K_e为渗流-管流耦合模型的等效渗透系数。

河西走廊疏勒河流域地下水资源合理开发利用调查评价

需注意的是，紊流条件下的等效渗透系数K_N与流速v及摩擦系数f有关，它是个随雷诺数而变化的量。另外，雷诺数Re的确定与流速v有关，而流速v又依赖于摩擦系数f，摩擦系数f的确定反过来又取决于雷诺数Re，因此必须采用迭代法来确定三者，进而确定等效渗透系数K_e。具体计算流程如图5-4。

二、混合观测孔水位的模拟

由于种种原因，几乎没有无混合观测孔的地区。什么是混合观测孔，混合观测孔中的水位（混合水位）如何形成，如何根据混合水位来确定初始水头的分布，如何利用混合水位求取各分层水文地质参数，等等，是“防止模拟失真，提高仿真性”（陈崇希，2003）的重要问题之一。然而MODFLOW（McDonald等，1988）对混合观测孔没有做任何分析、讨论，也没有混合观测孔的模块。简单地放弃混合观测孔的信息是不对的，因为“混合观测孔”只是孔口的流量为零，混合观测孔内却存在“抽水”与“注水”。

D.Sokol（1963）采用Thiem“影响半径”稳定井流模型证明混合观测孔的水位等于以导水系数为权重的各层水头的代数平均值。该论文存在如下几个主要问题：①Thiem“影响半径”模型是不能形成稳定流的（陈崇希，1966，1975，1983）；②Sokol隐含着假定，尽管两含水层的导水系数不等，流量也不同，但取其“影响半径”相等；③Sokol还隐含着假定，井管的阻力可以忽略不计，即井管处处水头相等。在这些假定条件下得出的结论是不可信的。

Hantush（1961）和Бочевер等（前苏联）（1961）在对承压含水层非完整井流的研究中，都认为观测孔中的水头降深s（r，l′，d′，t）反映该孔滤水管中各点降深s_p（r，z，t）的平均降深，即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：l′，d′是观测孔滤管顶点、底点的标高。本书称之为Hantush-Бочевер方程。

对于潜水三维流的研究，Neuman（1972）也认为，观测孔中的水头降深可视为滤管内各点降深的平均值，即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：r是观测孔至抽水井的距离；z₁、z₂是观测孔滤管顶点、底点的标高。

上述诸学者的见解与确定观测孔中水位的方法在水文地质界一直沿用至今，40多年来未见到异议。然而，上述观测孔中水位的确定方法，缺乏对形成机理的基本分析。简单分析之，三维流中的垂直观测孔，由于观测孔滤管中的水头不等，会发生垂向流动，而水的流动则会导致水头的变化，并引起井孔周围地下水的运动与水头再分布。因此，与测压计式观测孔不同，一般观测孔实质上是一部分滤管进水（抽水）而另一部分滤管出水（注水），即观测孔并不简单地反映含水层地下水的水头，而是它兼有“抽水”与“注水”作用的井孔，只是孔口的流量为零，井管内的“抽水”量与“注水”量绝对值相等而已（如果忽略微小的井筒储存量的变化）。观测孔非“观测”孔也！

显然，上述Hantush、Бочевер等以及Neuman提出的“积分平均水位”都未涉及水流的机理，是缺乏物理基础的纯数学方法。作者曾用“渗流-管流耦合模型”来模拟混合观测孔中的水位（混合水位）。算例的条件是：

一个均质、等厚、水平无限延伸的承压含水层，于顶面处有一微小的半球状井定流量抽水，抽水延续时间t=10.02d，在不同径距r处有不同孔径的完整观测孔（此即混合观测孔），模拟这些观测孔中的水头分布和流量变化。含水层的厚度M=100m，渗透系数K=100m/d，单位弹性给水度μ_s=0.00001m^-1，抽水流量Q=36000m³/d。模拟结果示于图5-5和图5-6。由图5-5可见，孔径≤0.1m的观测孔中，不同深度处的水头值明显存在差异，与Hantush-Бочевер的积分平均相差也很大。

图5-5 不同径距r观测孔中水头降深分布及Hantush-Бочевер公式计算值s_ga

图5-6 不同径距观测孔流量分布

由于抽水井（点汇）位于承压含水层顶面处，因此所有径距r处的顶部处的水头降深s大于底部，即底部的水头H大于顶部的水头，于是观测孔中的水从下而上流动。基于水流连续性原理，底部观测孔要吸水（抽水）而顶部观测孔要排水（注水），图5-6表示观测孔中的流量分布。图5-6a表明：径距21.19m，直径0.20m的观测孔中，大约在z=35m以下为“抽水”，35m以上为“注水”；观测孔中垂向最大流量超过500m³/d。这就是我们所说的，观测孔非“观测”孔也！

基于上述分析，本项研究采用“渗流-管流耦合模型”来模拟混合观测孔中的水位，使混合水位成为有用的信息，即用于初始水头分布和拟合求参。

三、泉流量动态的模拟

泉是地下水转化为地表水的主要形式之一，不同类型泉的出现及其流量动态，是水文地质工作者要重点分析的问题，因为这些数据为认识水文地质条件提供了极为重要的信息。泉是数值模型中主要的模拟要素；泉流量的预测是地下水资源评价、管理不可缺省的内容。

美国地质调查局（2001）建议采用MODFLOW软件的Drain模块迭代计算泉流量，实质上将泉作为第二类边界来刻画，即先不考虑泉的存在求解水头分布，再用泉所在格点的地下水水头与泉口标高之差乘以某比例系数（该系数缺乏物理意义）计算出泉流量，再以此流量置于该格点（第二类边界）重新求解水头分布，再次重新计算出新的泉流量，如此迭代直至收敛为止。当地下水水位低于泉口标高时，流量为零；当地下水水位高于泉口标高时，泉流量与水头差成正比，其比例系数必须通过流量的拟合来反算。

我国以往数值模型中关于泉流的刻画，采用的像抽水井一样输入（给定）流量。如此处理，预测怎么办?泉流量是未知的。近年来，我国也有类似上述MODFLOW建议的迭代算法。例如有的研究报告将鄂尔多斯大向斜的东侧山西省柳林泉（实为单斜自流斜地上升泉）刻画为平面二维流模型，进而用承压完整井的公式计算泉的流量，即Q=qs。这里存在几个问题：①柳林泉是上升泉，泉附近的钻孔愈深其地下水头也愈高，这是单斜自流斜地上升泉的特点，这种条件具有明显的三维流特征；②用承压完整井的公式计算泉的流量，特别是柳林泉的流量，不妥；③比例系数q在不同时间（不同条件）并不是常量。

“泉，必需将其放在水文地质体中去研究其成因类型，才能正确模拟仿真”。本模型提出的方法不同于MODFLOW等，是将泉口标高作为第一类边界条件（当泉流存在时），并将泉口格点与同层周边格点及下格点分别依达西定律（如果地下水流属于线性流，例如陈崇希等（1995），Chen Chongxi等（1996），陈崇希等（2003）；或采用非达西定律，例如陈崇希（1995），成建梅等（1998），ChengJianmei等（1998））和水均衡原理建立关系，整个计算方法建立在流动机理之上，因此由模型运转直接产生泉流量，无需通过迭代求解（且有的情况下迭代并不收敛）。正因为此，在模型识别阶段才有可能将泉流量作为拟合对象。这一点十分重要。

具体地说（图5-7），根据达西定律，泉流量Q为水平向的流量与垂向上的流量之和，如果以水平向的流量为主，则为下降泉，否则为上升泉。具体可通过下式计算：

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：Q为泉流量；K_zsp为泉流区域的垂直渗透系数；H_b为泉口下层结点的水头值；Z_sp为泉出露标高；D_zsp为泉流过程中的渗流长度（含水层厚度）；S为泉口结点控制面积；C_A为泉颈修正系数；C_s为泉径修正系数；T_ij，T_ik为流段ij和ik的平均导水系数；，，，为线段的长度；h_i，h_j，h_k为结点i，j，k的水头值；e为结点i周围的均衡单元，如ipbq四边形。

图5-7 泉的流量计算示意图

在以往的地下水资源评价数值模拟中，没有将泉的流量作为拟合要素。实际上，象MODFLOW等的做法也不可能将泉的流量作为拟合要素，因为他们在模型识别其间要求解一个缺乏物理意义的系数。

本模型选择了3个典型点泉和22条线泉（泉沟），模拟其流量动态，并将其实测流量作为一个重要的拟合对象用以求取含水系统的水文地质参数和预测泉流量的动态。如此将会提高模型的仿真性。

四、大气降雨、地表水等入渗补给潜水滞后性的刻画

降水、地表水、渠系等入渗补给潜水是地下水的重要补给来源，若处理不当，则会导致模型失真。因此，正确的刻画入渗补给的滞后性对于模型仿真性的提高具有重要的意义。

目前地下水资源评价中降水、地表水、渠系等入渗补给量的计算主要采用下列两种方法：

1）未考虑滞后的入渗系数法。该法对于小埋深区域大体上可以用来刻画入渗对潜水的补给，但有些地区潜水埋深较大，若不考虑降雨、地表水等入渗补给的滞后特征，则会导致模型失真。

2）平移滞后入渗系数法。有的研究报告把潜水位出现高峰时间与降雨高峰出现时间的时间差，作为统一的滞后时间。例如若该时间差为3个月，则1月份的降雨量都在4月份入渗补给，2月份的降雨量则在5月份入渗补给；若3月份无降雨，则6月份无降雨入渗补给。如此处理滞后补给也存在明显的不合理性。因为潜水位高峰出现的时间，是所有影响地下水位动态诸因素作用的总和，除降雨补给外，还有地下水蒸发，灌溉入渗，越流和地下水开采等因素作用下产生地下水不平衡流动的结果，特别是地下水开采这个人为因素在本区已经成为地下水动态的主要控制因素，它已经完全掩盖了地下水的天然动态，因此用该法研究降雨入渗补给的滞后性是不恰当的。即使这一地区完全处于未开采状态，甚至忽略除降雨以外的其他因素对地下水动态的影响，仅就降雨入渗补给单一因素对潜水位动态的影响而言，该法也存在明显的问题。降雨过程是间断的，不连续的，但它对潜水的补给则是连续的，这是因为地下水流动，特别是非饱和流动，其介质具明显的储存性而有调节功能，从而使补给具有滞后性及连续性。一年内虽然降雨仅有几十天，其余时间均无降雨，但在潜水位深埋区，补给却天天作用着（浅埋区由于存在强蒸发作用而使情况复杂化了），这种连续补给作用是不可能用一个时间差的平移来刻画其滞后性的。

上述“入渗系数法”在潜水小埋深区有一定的适用性。但由于本区部分地段潜水埋深较大，如昌马洪积扇顶部潜水水位最大埋深达290m上下。一个月的降雨量乘以入渗系数的水量不可能在当月内全部补给其下的潜水含水层。因此，不考虑降雨补给的滞后特征，则会导致模型的失真。

本研究采用“降雨补给滞后权系数法”（陈崇希等，1991，1998）来刻画降雨入渗补给。该方法既能基本反映客观实际情况，又很实用。

根据地下水非饱和流动理论，一次降雨过程对潜水补给量的分布曲线如图5-8所示。它是一单峰曲线。曲线的形态取决于潜水位的埋深和包气带的岩性，当然还与包气带的初始含水率的分布有关。我们将后面的因素固定起来，专门讨论前面因素的作用。

图5-8 一次降雨对补给量的分布曲线图

我们以包气带地下水非饱和流动理论为基础，考察不同水位埋深和岩性条件下，某时段（例如月、旬等）降雨对该时段及其后各时段补给强度的分布特征。为了使问题的讨论能与地下水资源评价数值法相匹配，将连续的特征曲线改为离散点的形式（因为数值法在时间上是离散的）。不失一般性，我们将时段取为1个月。

对于潜水位埋深较大和（或）包气带渗透系数较小的条件，入渗强度分布是一单峰函数，当潜水位埋深很小时，可能是单调减函数（图5-9）。若在0月有某降雨量RA（mm）作用，则该月及后各月降雨入渗补给量之和∑RE_k（mm）就是该时段降雨量形成的全部补给量。该值∑RE_k与0月的降雨量RA之比就是（总）入渗补给系数α。各月份的入渗量与0月降雨量RA之比，则是0月降雨在各月的入渗补给系数α_k。即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

图5-9 降雨对不同埋深（岩性）条件的入渗补给历时分布曲线图

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：RA为0月的降雨量，mm；RE_k为0月的降雨在k月的入渗补给量，mm；RE为（总）补给量，mm；α_k为0月降雨在k月的入渗补给系数；α为（总）入渗补给系数。

如果令图5-9上0月的降雨量为100mm/月，则该图纵坐标转变为以百分数表示的月入渗系数α_k（k=0，1，2，3，…）。

月入渗补给系数组已能很好地反映降雨入渗补给滞后性的特征，为了使刻画滞后性的特征数组更通用，我们定义：

河西走廊疏勒河流域地下水资源合理开发利用调查评价

因此

河西走廊疏勒河流域地下水资源合理开发利用调查评价

该式表明ω_k满足权的定义，所以我们称ω_k（k=0，1，2，…）为滞后补给权系数（组）。

引入滞后补给权系数（组）的优点是，不同入渗系数的诸分区（不同区，同月的入渗补给量不等），有可能取用同一滞后补给权系数来计算各月的补给量，也就是说，滞后补给权系数更具通用性。于是我们的任务转到寻找一个能满足上述条件的特征函数，即：

1）它们是单峰函数或单调减函数，并随峰值（或最大值）出现时间的后延，其峰值（或最大值）减小；

2）各权系数之和为1；

3）适应性较强，能较好地与不同潜水位埋深，不同岩性条件下滞后入渗补给量相拟合。

经反复分析研究，我们确定采用下面的离散函数来表征滞后入渗补给权系数的分布。

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：ω（j-i）为i月降雨，j月补给的权系数；h为潜水位埋深，m；B为与包气带岩性等因素有关的经验系数，m；可称为埋深区间系数。该系数是包气带渗透系数的增函数。

该式中h是已知的，仅有一个与岩性等因素有关的经验系数B是待求的。这个系数可以通过模型识别来确定，这样就可以利用上式计算出滞后入渗补给权系数。各区（总）入渗系数也是通过模型识别来确定的，由此可计算出在某月i降雨量RA（i）作用下，该月及其后各月j的补给量RE（j-i），即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

至此，我们分析了某月降雨量对该月及其后各月入渗补给的贡献。对于实际问题——多月降雨对某月j的补给量问题，为了简化计算，我们近似用该月及其以前各月i降雨量对潜水补给贡献的代数和来计算该月j获得的总入渗补给量，即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：RES（j）为j月及其以前各月的降雨量在j月的总入渗补给量，mm；RE（j-i）为j月降雨对j（j＞i）月的入渗补给量，mm；α为入渗补给系数；ω（j＞i）为i月降雨对j（j＞i）月的滞后补给权系数；RA（i）为i月的降雨量，mm。

依式（5-17），获得的降雨入渗补给滞后权系数曲线共75条，它对应于75个潜水位埋深区间，每个区间埋深差Bm（即埋深区间系数）即每增加埋深Bm，换一条权系数曲线的序号）。图5-10中仅绘出12条典型的权系数曲线，而各两条曲线之间，尚有5～7条曲线图中并未绘出。

降雨补给滞后权系数法的提出和使用，不仅大大提高了模拟的仿真性，而且使用十分方便。地表水入渗补给滞后性的机理与降雨入渗的类似，也采用上述方法刻画。

图5-10 降雨入渗补给滞后权系数曲线

五、参数的分区

本模型含水介质厚度大，垂向上共分6个模拟层，即便是同一岩性（如亚粘土）而处于不同埋藏深度，由于其固结程度不同而导致其水文地质参数不相等。这种情况下如何给定初始参数估计值的空间分布又是一个难题。若按常规惯用的方法，则同一岩性根据不同层位（不同埋藏深度）要给出不同参数分区，如此的参数分区数过多，给调参带来很大的工作量。

对此问题，本模型采用如下方法，处于不同埋藏深度同一岩性的K_h值作为同一分区，但其参数值随埋深的负指数衰减。这样，每个参数分区只要1个零埋深的水平渗透系数K_h0和1个衰减系数γ刻画，即K_h=K_h0×e^-γh（h为埋深值）。其他参数，如K_z值、μ_d值和μ_s值的空间分布也按同样方法给定。如此处理，就大大地减少了模型参数的数量，相应地也减少了模型调参的工作量。

六、自流井的模拟

自流井普遍地存在于我国西北内陆盆地，冲洪积扇的前缘洼地，大型岩溶上升泉的周围以及平原区的深部承压含水层等区段。我国至今还存在大面积的地下水自流区（尽管由于大量开采地下水而缩小了其自流面积）。因此模型中如何刻画自流井是一个重要的水文地质问题。然而，美国地调局推出的MODFLOW软件尚无模拟自流井流量的模块，这似乎不是疏忽，而是存在一定的难度。常规的地下水流模型已经解决不了这类问题。至今，对于数值模型中的自流井，仍是无可奈何地采用人为地给定其自流量。然而，自流井与抽水井不同，后者人们一般可以人为给定抽水流量，而井中水位是作为其响应；当井中水位降至井底时，人们就不能自由地抽取了（尽管如此，大多含有抽水井的数值模型，未能客观地将“井中水位作为抽水量的约束”来认识）。然而，自流井就不能人为地给定自流量了，因为自流量是自流井对周围环境（气象、水文、地下水的开采等因素）的响应，自流井的流量是属于预测范畴，是输出信息，而不是输入信息。

本项目采用“渗流-管流耦合模型”模拟自流井，只要在自流井的井口处设置水头已知边界，即水头等于井口标高，即可模拟出自流井的流量。该法已在“陕西渭北东部岩溶水资源数值模拟研究”中获得成功的应用（陈崇希等，2003）。

七、初始水头的形成

初始水头的分布是地下水不稳定流数值模拟不可缺少的条件。通常利用一定数量观测孔的水位通过插值获得各结点（格点）的初始水头值。然而，当研究区范围较大，观测孔较少，观测孔所在的层位不一，而且含有混合孔时，难以用插值的办法获得每一模拟层的初始水头分布。

本项目采用“参数-初始水头迭代法”（陈崇希、林敏等，1991；陈崇希、裴顺平，2001），来确定各层的初始水头分布。具体处理办法是：取模拟的初始时间为2003年7月，用水文地质参数初始估计值从2001年7月开始模拟，至2003年7月将模拟所得的该月水头值与实测值拟合，粗略确定模拟时段的初始水头（2003年7月）；再以此初始水头按常规方法求参，如此迭代求解。当拟合达到一定程度之后，必须全面地检验模型是否符合基本规律，否则要再次重新调整，直至满足要求。上述方法虽然工作量很大，但效果良好。

‘陆’ FEFLOW地下水模拟软件系统简介

FEFLOW（Finite Element subsurface flow system）是由德国水资源规划与系统研究所（WASY）历时二十多年的研究，开发出来的地下水流动及物质迁移模拟软件系统^［9］。该软件提供图形人机对话功能、具备地理信息系统数据接口、能够自动产生空间各种有限单元网、具有空间参数区域化、快速精确的数值算法和先进的图形视觉化技术等特点^［22］。

FEFLOW最初是用FORTRAN Ⅳ 编写的，并且作为典型的批处理程序用于大型计算机，在20世纪80年代，FEFLOW推出了由C语言源程序改写的第二个版本^［21］。此后，在理论研究和实际问题的处理上，经过了不断的发展、修改、扩充、提高，日趋完善。从20世纪70年代末至今，FEFLOW 经过了大量的测试和检验，成功地解决了一系列与地下水有关的实质性问题，如判断污染物迁移途径、追溯污染物的来源、海水入侵等等，是迄今为止世界上功能最齐全、技术最为先进的三维地下水模拟分析软件。目前，最新的FEFLOW版本是VERSION5.3，由C语言和C++语言编写而成。

5.3.1.1 FEFLOW 的应用领域

FEFLOW的应用领域包括^{［22、24］}：

模拟地下水区域流场及其他地下水资源规划和管理方案；

模拟矿区露天开采或地下水开采对区域地下水的影响及其最优对策方案；

模拟由于近海岸区抽取地下水或者矿区抽排地下水引起的海水或深部盐水入侵问题；

模拟非饱和带以及饱和带地下水流及其温度分布问题；

模拟污染物在地下水中的迁移过程及其时间空间上的分布规律（可用于分析和评价工业污染物及城市废物堆放对地下水资源和生态环境的影响，研究最优治理方案和对策）；

结合降雨—径流模型进行模拟“降雨—径流—地下水”的系统问题（可用于研究水资源合理利用、管理以及生态环境保护方案等）。

5.3.1.2 系统输入特点（建立模型）

通过标准数据输入接口用户既能直接利用已有的GIS空间多边形数据生成有限单元网格，也可以用鼠标设计和调整网格几何形状，增加和放疏网格密度^［22］。在建立水流场和迁移模型时，用户不仅能够视具体情况定义第一、第二和第三类边界，而且可以对边界条件增加特定的限制条件，以避免非现实的数值解。用户也能够直接定义多含水层中的抽水和注水井边界条件。所有边界及附加条件既可设置为常数，也能定义为随时间变化的函数。已知的边界及模型参数可以按点、线或面的形式直接输入。对离散的空间抽样数据进行内插或外推（数据区域化），FEFLOW提供克里格法（Kriging），阿基玛（Akima）和距离反比加权法（IDW）。输入数据格式既可以是ASCII码文件，也可以是GIS地理信息系统文件。FEFLOW支持ARC/INFO点、线、面的广义数据格式，ArcView形状数据格式，DXF格式，Tiff图形以及HPGL数据格式。

5.3.1.3 系统模型求解特点

FEFLOW具有齐全的地下水模拟功能^［22］：

三维空间模型，二维平面，二维剖面或者轴对称二维模型；

非稳定流或稳定流模拟；

多层自由表面含水系模拟，包括滞水（perched water）模拟；

化学物质迁移及热传递模拟，包括温度盐分（thermohaline）迁移模拟；

可变密度流场模拟（盐水或海水入侵问题）；

非饱和带流场及物质迁移模拟。

FEFLOW采用加辽金法为基础有限单元法，并配备若干先进的数值求解法来控制和优化求解过程：

快速直接求解法，如 PCG，BICGSTAB，CGS，GMRES 以及带预处理的再启动ORTHOMIN法；

灵活多变的up-wind 技术，如流线up-wind，奇值捕捉法（Shock capturing）以减少数值弥散；

皮卡和牛顿迭代法求解非线性流场问题，自动调节模拟时间步长；

模拟污染物迁移过程包括对流，水动力弥散，线性及非线性吸附，一阶化学非平衡反应；

为非饱和带模拟提供了多种参数模型如指数式，Van Genuchten 式和多种形式的 Rich⁃ard 方程；

垂向滑动网格（BASD）技术处理自由表面含水系以及非饱和带模拟问题；

适应流场变化强弱的有限单元自动加密放疏技术，以获得最佳数值解；

实时图形显示模拟非稳定流过程中观测点水头和污染物浓度的动态变化值；

非稳定流模拟计算可以随时暂停，以便用户显示和分析中间模拟结果；

开放性外部程序接口，以便用户在FEFLOW 系统中连接和使用自己的程序模块。

5.3.1.4 系统结果输出及显示

FEFLOW的计算结果既有水位，污染物浓度及温度等标量数据，也包括流速，流线和流径线等向量数据。模型参数和计算结果既能按ASCII码文件，GIS地理信息系统文件，DXF或HPGL文件输出，又能在FEFLOW系统中直接显示和成图。FEFLOW提供了其他任何地下水模拟软件都无法比拟的、丰富实用的图形显示和数据结果分析工具^［22］。

5.3.1.4.1 先进的图形视觉化及数据分析技术

有限单元网，边界条件和模型参数的三维可视化及显示；

标量数据的三维彩色（透明或灰度）等势面显示及其二维平面彩色或等值线显示；

三维地下水流径追踪，流动时间及流速动画显示（包括其二维平、剖面投影或二维平面追踪）；

三维体截段的空间显示和三维交叉剖面组的空间显示；

三维图形的任意旋转，二维、三维图形的放大或缩小；

总体和局部水量平衡分析（包括任意几何多边形内的水流通量分析）；

计算和图形显示通过各种边界的水通量，物质通量及其在特定时间区间内的积分量；

借助FEFLOW 提供的XPLOT 或FEPLOT 程序模块可以直接设计和打印各种成果图件。

5.3.1.4.2 多功能的先进技术

能处理多种点阵式和向量式背景底图；

充分利用已有的ARC/INFO GIS地理信息系统数据产生有限单元网，设置边界条件和参数；

自动生成各种有限单元网并分析其几何特性；

模型空间分布参数的内插和区域化；

对模型边界条件及参数设置随时间变化的附加限制及规定；

实时显示地下水非稳定流场、温度场及污染物迁移模拟结果；

自动调整有限元网格密度，以获得随流场时空变化的最佳数值解。

5.3.1.5 应用方法

利用FEFLOW进行地下水模拟时，主要有以下几个步骤^［9］：

（1）确定模型的基础结构：为了定义有限元模型的外边界，首先添加背景图，用之作为基础，产生一个数字化的超级单元网格，这个超级单元网格将作为模型的基本结构；通过进入Edit菜单下的Mesh Generator menu子菜单，利用已生成的超大元网格可以自动生成有限单元网格，网格的数目可以自己指定，可以方便地调整网格的几何形状，增加和放疏网格大小等。

（2）确定模型的三维结构：根据模型的需要，可以进入Dimension 菜单，来定义片（Slices）和层（Iayers）的数目以及每个层的厚度，若为二维模型则此步忽略。

（3）输入模型的初始条件、边界条件及参数：进入Edit菜单下的Problem Editor Menu子菜单定义模型的参数。从问题编辑器的下一级子菜单Flow Data进入Flow initial，通过激活database菜单导入 FEFLOW 支持的 Arc/Info、数据格式、ArcView 形状数据格式、DXF格式、Tiff图形以及HPGL数据格式等，利用模型的自动插值计算功能进行自动插值计算，得到模型区域的流场等值线图；在Flow boundary菜单下定义模型的边界类型，并通过数据库或直接键入值来输入属性值；在Flow materials菜单下编辑所有模拟地下水流问题所需要定义的参数，可以通过导入数据库，或通过单元或面状分区在模型中直接键入参数值。对每个层来说，模型的初始条件、边界条件和水文地质参数可以分别输入，如果每层的条件相同，可以利用层间的COPY功能来实现赋值。

（4）模型模拟：由模拟程序外菜单Run进入模拟程序Simulator，敲击Run simulator按钮，开始模拟。FEFLOW可以自动生成并详细列出有关井的水头值的信息窗口，中间可以通过按住鼠标右键不放来停止模拟的进行，重新模拟用（Re）Run Simulator开始。

（5）分析模型结果和数据输出：模拟程序后处理菜单Postprocessor允许用户赋值、分析和输出计算结果。利用FEFLOW后处理菜单强大的数据分析功能，可以实现有限单元网、边界条件和模型参数的三维可视化及显示；标量数据的三维彩色透明或灰度等势面显示及其二维平面彩色或等值线显示；三维地下水流径追踪，流动时间及流速动画显示包括其二维平面、剖面投影或二维平面追踪；三维体截段的空间显示和三维交叉剖面组的空间显示；三维图形的任意旋转，二维、三维图形的放大或缩小；总体和局部水量平衡分析（包括任意几何多边形内的水流通量分析）；计算和图形显示通过各种边界的水通量，物质通量及其在特定时间区间内的积分量；借助FEFLOW提供的XPLOT或FEPLOT程序模块可以直接设计和打印各种成果图件。

5.3.1.6 应用FEFLOW 时需要注意的问题

应用FEFLOW时需要注意以下问题^［9］：

（1）在定义模型的三维结构时，先允许用户通过设置顶片的虚拟高程和层厚度来分配层，例如，已通过Layer configurator设置顶片高程为1000m，其下方的其他层间距皆为100m，然后可以通过设置相关的高程值将每片拉到其实际的位置。这时，为了避免交叉，应先设置最底下一层的高程数据。

（2）使用FEFLOW时，在各含水层初始水位的导入、地层标高数据的导入插值时特别需要注意一些关键的问题，如存在不合理的插值数据，输入模型中可能会导致模型不能运行，或是出现不合理的运算。

（3）FEFLOW菜单中的参数子菜单比较多，参数的名称、单位和意义在进行地下水流数值模拟时需要特别注意，只有了解其在地下水流计算中的意义，然后相对应给定参数值，计算结果才会准确无误。

（4）进行模型边界属性赋值时，注意边界上流入为负，流出为正。

‘柒’ 埕北古4井区中深层储集层预测及开发对策

宋美虹季雅新王玉芹杜玉山

摘要针对埕北古4井区东营组油藏特有的地质特征，在深入分析埕岛油田东营组相似区块开发经验的基础上，应用储集层预测技术、数值模拟方法以及初步的经济评价，对该井区的储集层分布进行了预测；对其关键开发技术政策进行了优化研究；并以此为依据进行了开发方案部署。

关键词埕北古4井区储集层预测相似区块类比数值模拟经济评价开发方案

一、引言

目前，海上产能建设的主阵地已由浅层转移到中深层。中深层油藏地质特征复杂，地震资料反射能量弱，中高频率损失严重，信噪比与分辨率都较低。适合于浅层的储集层预测技术在中深层已不适用。因此，发展完善中深储集层地震预测技术，对落实其石油地质储量，提高总体开发效果和经济效益，具有重要指导意义。

埕北古4井区东营组是主力含油层系。油层埋深在2900m以下，属于中深储集层。为了描述储集层、落实石油地质储量、进一步指导开发方案的编制、降低开发井的部署风险，通过多方调研，查阅了国内外有关信息资料，对本区先后采用TRP软件和JASON软件进行了地震储集层预测研究。实践证明，这两次攻关研究是较成功的。

由于中深层油藏的油藏类型、储集层特征及流体性质与主体已开发的浅层馆上段油藏不同，浅层的开发技术政策不适用于中深层油藏。在对埕岛油田东营组已开发区进行深入分析的基础上，针对本区含油井段较短、油层少且主力层突出、油水关系简单、储量规模不大等地质特点，制定了开发原则。通过数值模拟和经济评价，对本井区的开发方式、布井方式、采油速度以及水平井段长度等进行了研究，确定了以天然能量开采为主、辅助于注水补充能量开发的开发方式，优选出定向井与水平井组合的布井方式。方案设计总井数5口，其中定向井3口，水平井2口。动用储量556×10⁴t，预计建成年产能力18×10⁴t，开发15年，可累积生产原油102×10⁴t。

对埕北古4井区中深层特有的地质和油藏特征、储集层预测及开发对策做了深入研究，总结出了一套相关的方法。该研究成果对具有类似地质特征的油田新区储集层预测及方案设计具有参考价值。

二、工区概况及其油藏地质特征

1.工区地理位置及勘探现状

埕北古4井区位于埕岛油田东北部，西邻胜海古2及胜海古3井，南与埕北8井相接，水深15～20m，构造位置位于埕宁隆起埕北低凸起东斜坡下第三系超覆带。

1996年1月10日完井的胜海8井为该井区第一口探井，完钻井深3600m，完钻层位中生界，电测解释油层2层26.7m，均为东营组。经测试，在2021.3～3052.0m井段获日产油224t，气22453m³。此后，又相继完钻3口探井。埕北古4井区4口井6个层段试油，3口井四层段获工业油流。埕北古4井区有3口井试采。

2.地层层序及含油层系

该井区自下而上钻遇的地层有古生界、中生界、下第三系沙河街组、东营组、上第三系馆陶组、明化镇组及第四系平原组。发现了古生界、中生界、沙河街及东营组四套含油层系。其中，东营组是埕北古4井—胜海801井区的主力含油层系。

3.东营组沉积特征

东营组其下部为湖相沉积，中间为扇三角洲前缘亚相沉积，上部为扇三角洲平原亚相沉积。其底部发育大段泥岩夹薄层砂岩，中间发育大段厚层砂岩，上部为砂泥薄互层，具有“底超顶剥”的特点。埕北古4井—胜海801井区东营组构造位置比较低，地层发育相对比较全，只有埕北古4井下部缺失东营组部分地层，其他3口井底部地层均发育齐全。埕北古4井、胜海801井及胜海8井顶部有少量地层被剥蚀。区内东营组地层厚度为610～890m。其中，埕北古4井和胜海801井厚度小，胜海8井及胜海10井厚度大。

4.储集层特征

埕北古4井—胜海801井区东营组可分为8个砂层组，油层主要分布在6、7砂层组，油层埋藏深度为2914～3430m，含油井段长516m，平均单井钻遇东营组油层14.9m/3层。最大单井有效厚度29.5m，最小4.0m。其中，埕北古4井钻遇油层最多，共29.5m/7层，均为一类有效厚度；胜海801井钻遇油层最少，共4.0m/1层，为一类有效厚度；区内最大单层有效厚度14.5m，最小单井有效厚度1.1m。油层发育主要受岩性控制，其次受断层控制。只有与断层相接触的砂体才可能形成有效圈闭而含油。

三、储集层预测研究

1.原始地震资料品质分析

本次储集层预测处理地震资料面积约60km²。涉及5口井（胜海10、埕北古4、胜海8、胜海801、埕北古403）。所用的地震资料时窗为1500～3500ms，采样间隔2ms，三维网格为25m×25m。

地震叠偏数据体的分辨率、信噪比、保真度等品质分析如下。

（1）分辨率

目的层的平均速度取3000m/s，可分辨厚度为λ/4；目的层顶部视频率约30Hz，分辨厚度约25m；目的层中部视频率约26Hz，分辨厚度约29m；目的层底部视频率约22Hz，分辨厚度约36m；

（2）信噪比

总体看，该区地震资料信噪比较好，尖灭点、断点、超覆沉积现象比较清晰，但不足之处是剖面偏移划弧现象严重，造成某些断点不清和偏移干扰等负效应。

（3）保真度

经过偏移处理的地震资料，数据格式是32位浮点，2ms采样，数据体能量较均衡，资料有一定的保真度^[1]。

2.储集层预测采用的方法

为了描述储集层、落实石油地质储量、进一步指导开发方案的编制、降低开发井的部署风险，先后采用两种方法进行了地震储集层预测研究^[2]。先于1998年底，采用TRP软件；后于1999年3～5月，采用JASON软件。

1)TRP软件

（1）反演原理

由测井资料给定的初始阻抗模型Z(t）与从地震资料提取的实际子波W(t）正演模拟得到当前道的合成记录。在地质模型、地震特征约束下，通过合成记录与实际记录的相关对比，经过反复迭代来调整当前道的波阻抗模型。当两者误差满足要求时，对应的阻抗模型即为当前道的波阻抗反演结果。

（2）处理流程

TRP软件为井约束下的高精度三维储集层参数反演技术。其处理过程主要包括测井资料处理、建立单井波阻抗模型、建立精细地震地质模型、多井约束三维波阻抗反演等。

测井资料处理 对各井测井资料进行环境校正，消除井径、泥浆滤液及压实作用等对测井曲线的影响。经标准化处理后，井之间波阻抗的相对差异则有较大程度的削弱，且横向的非均质性也得到较好的保持。

建立单井波阻抗模型 建立较准确的井中时深关系及井中波阻抗与岩性的对应关系。利用井的速度和密度曲线得到井的波阻抗，进而得到深度反射系数；根据井的时深关系求得时间域的反射系数；然后从井旁道提取子波，子波与反射系数褶积得到井的合成记录；依据合成记录与井旁道相关系数最大原则，扫描出最佳井旁道。同时，也求得了最佳时深关系。再通过迭代调整得到最佳井旁道的波阻抗作为井的波阻抗，最佳井旁道就是井的对应地震道——种子道，也就是用遗传算法进行波阻抗外推的母体道。绘出各井的波阻抗图进行对比分析，如发现不合理之处，返回上一步调整处理，直到合理为止。

三维地震地质模型解释和处理 三维地震地质模型处理包括层位处理和断层处理。首先从井的标定层出发，对该区的层位进行详细的追踪解释，解释出层位和断层；其次，通过插值得到断面数据；对层位加断层边界进行内插，得到该层的层位平面图。对该区所有层都进行如此处理，结果即为三维地震地质模型。

三维波阻抗反演 三维波阻抗反演首先要形成控制文件，然后进行单井波阻抗反演，最后进行临接区反演。

控制文件确定了一口井的反演范围和顺序，一口井控制的范围是一个多边形区域，一般以断层为边界，以目的层为趋势，在井间留出20道以上的过渡区。

根据单井控制文件，从井对应的种子道出发，一圈一圈地向外反演。要反演的地震道从以自身为中心的小面元中选取最佳初始波阻抗模型，在地质模型、地震特征约束下，用迭代法调整所选定的初始波阻抗模型，当波阻抗对应的合成记录与当前地震道的相关系数达到80%以上时，认为得到了当前道的波阻抗结果。因为考虑了地震波场在各个方向的分布与变化因素，反演结果稳定可靠。

临接区包括断层区和井间过渡区。对于临接区的道，首先从八个方向在要求的范围内查找已反演的道，在查找过程中，如果那个方向遇到断层，或者在要求的范围内找不到已反演好的道，则忽略这个方向；然后用这些道插值出临接区道的初始波阻抗，再通过迭代得到该道的波阻抗结果。

（3）软件反演特点

按传统算法以井点标准波阻抗作为种子道分井区反演，利用完钻井的资料作为约束条件，提高了井周围预测精度；储集层成层性显示较好，储集层边界、断层清楚；反演剖面分辨率较高。

2)JASON软件

（1）反演原理及处理流程

本次反演工作采用了JASON软件中的稀疏脉冲反演方法。其基本假设是反射系数是稀疏的。该方法适用于区内井数较少的开发准备阶段。其主要优点是能获得宽频带的反射系数，从而使反演得到的波阻抗模型更趋于真实。

稀疏脉冲反演的主要过程是：通过最大似然反褶积求得一个具有稀疏特性的反射系数序列；通过最大似然反演导出宽带波阻抗。

（2)JASON反演的特点

子波的选取 在地质框架模型基础上，利用多井估算多个子波，由控制点估计的子波进行内插得到空变子波。能够监控子波的波形、相位和频谱，也可以监控各种各样的子波所产生的效果。使用空变子波合成的地震记录与地震资料的相关性好，从而达到最优效果。

宽频带的地震资料反演 JASON反演所导出的结果是一个宽频带的反射系数序列和宽频带的波阻抗数据，其低频分量是在地质框架模型基础上利用所有层速度建立低频模型，并与反演结果道道相并而得到的。从而保留了测井曲线的主要地质特征。

根据工区内所有完钻井建立整体三维模型，进行整体反演。

3.反演结果和砂体预测

根据区内已完钻的4口井分析，东营组5砂层组以上储集层很发育，但基本不含油，主力油砂体分布于6、7砂组。东营组的储集层主要沿斜坡带分布，高差大，空间速度变化大。依据本地区钻井资料及地质特征，结合两套软件处理的方法原理，在解释过程中，对剖面的色标不断加以微调，尽量准确反映主力油砂体，从而达到最佳预测效果。在反演剖面上，不同深度段的储集层颜色不完全一样。4口井储集层段的速度范围是3300～4100m/s。

结合储集层沉积模式，对两套波阻抗反演处理成果分别进行了地震储集层精细解释，重点描述埕北古4井区砂体。

TRP方法共解释了4个砂体，其中有井钻遇油砂体2个（7²、7^3～5），预测新砂体2个。

JASON方法共解释埕北古4井区3个砂体。在JASON软件处理的剖面上，埕北古4井7²及7^3～5中间的隔层反映较差，将7²、7^3～5作为一个砂体进行解释。另外，预测出了对应于TRP软件解释结果的两个新砂体。

4.预测效果分析

两銮软件预测的各砂体平面展布形态及面积大致相近，JASON软件预测各砂体的面积、厚度比TRP软件预测的要大些；两套软件反演结果与井的吻合都较好。从过井剖面上看，TRP反演结果的分辨率比JASON略高；本次反演所利用的井较少，在一定程度上影响了离井较远地区反演结果的精确性，在将结果应用于整个工区时还应该结合其他的资料进行综合分析，以提高决策的精确性。

四、关键开发技术政策研究

1.相似区块开发效果初步分析

到1999年9月底，埕岛油田东营组共上报Ⅲ类探明含油面积8.1km²，石油地质储量1429×10⁴t，主要分布在埕北11、12、21、斜101、35、151等6个区块，已于1993—1994年全部投入开发，累计建成产能30×10⁴t，累积产油132.1×10⁴t，采出程度9.3%。

埕岛油田东营组各区块的地质特征差异较大。针对本区地质开采特征，首先对东营组已开发区块进行筛选，即从含油层位、沉积类型、储集层物性及原油流体性质等方面进行对比分析，埕北古4井区与已开发的埕北21、151块比较相似，为了总结已开发区块经验，指导埕北古4井区东营组开发部署，对相似区块的关键技术政策进行系统分析。

（1）天然能量开发效果

埕岛油田东营组相似区块油藏地饱压差大（17.3MPa），具有较活跃的边底水，弹性产率为182×10⁴t/MPa。根据石油天然气行业标准，当采出程度1%时地层压降小于0.2MPa，弹性产量比大于30，为天然能量充足。而东营组相似区块采出程度1%时的地层压降为0.05MPa，弹性产量比为111，天然能量充足。

至1999年11月底，埕北21、151块上报探明储量385×10⁴t，共投产4口井，单井日产油能力29t，综合含水65.5%，年产油3.6×10⁴t，累积产油37.6375×10⁴t，采出程度9.8%，弹性采出程度较高。

（2）初期单井产量高，产量递减较快

埕北151块、埕北21块试油期间采油指数4.3t/(d·MPa·m），投产初期单井产量较高，平均日产油能力为113t，采油指数4.0t/(d·MPa·m）。

东营组投产井均依靠天然能量开发，产量递减较快。埕北21块由于只有埕北21井一口井生产，单井控制储量大，自1994年2月投产以来，产量一直处于相对稳定阶段，1998年5月进入产量递减阶段，递减较快，年递减率为27.6%；埕北151块产量递减较快，年递减率为32.4%。

（3）厚层块状油藏开采效果明显好于多层层状油藏

埕北151块油层层数多，平均单井钻遇油层4层18.7m，单层厚度小，油层连通性好，油水关系复杂，属典型的层状油藏；埕北21块油层单一，单层厚度大，钻遇油层1层31m，有效厚度29m，油水关系简单，属厚层块状油藏。

从相似区块的开采效果来看，属于块状油藏的埕北21块明显好于属于层状油藏的埕北151块。截止1999年9月底，埕北21井日产油63.1t，综合含水54.5%，累积产油26.6×10⁴t，采出程度11.3%；埕北151块开井3口，日产油77t，综合含水57%，区块累计产油10.62×10⁴t，采出程度7.1%。

（4）单井控制储量可以适当放大

埕岛油田东营组已开发区单井控制储量50×10⁴～235×10⁴t，平均65×10⁴t，略大于馆陶组单井经济极限控制储量（62×10⁴t），小于主体北开发区馆上段实际单井控制储量（81×10⁴t）。根据东营组各区块的开发效果，单井控制储量较大的埕北21井区（235×10⁴t）开发效果明显好于其他区块。

东营组较馆陶组储集层物性差，但原油性质较好，流动系数1487×10^-3μm²·m/(mPa.s），大于馆陶组（551×10^-3μm²·m/(mPa·s）；东营组油藏埋藏较深，岩石压缩性小，压力传导较快，导压系数5.56μm²·MPa/(mPa.s）大于馆陶组1.42μm²·MPa/(mPa·s）；东营组油藏埋藏较馆陶组深，单井钻井成本较大，因此，东营组油藏单井控制储量可以适当放大一些。

2.关键开发技术政策研究

为了合理地制定埕北古4井区东营组油藏的开发技术政策，在借鉴相似区块的开发经验的基础上，对本区又进行了数值模拟研究^[3]。

1）三维地质模型的建立

数值模拟目的层为埕北古4井实际钻遇的7²和7³5两个油砂体，模型区叠合含油面积3.58km²，石油地质储量376×10⁴t。

根据实钻井资料及储集层预测结果，建立了每个砂体的顶部深度、砂层厚度、有效厚度、渗透率、孔隙度等参数场。

在对相似区块和本区现有资料深入分析的基础上，结合相关图版及经验公式[4,5]确定了本区数值模拟所需的油藏工程参数，建立了岩石、流体等模型。

2）数值模拟方案设计

埕北古4井区油层较少且主力层突出、油水关系简单、储量规模不大。根据数值模拟研究内容（开发方式、不同布井方式、采油速度以及水平井段长度等）进行了数值模拟方案设计。首先，在叠合有效厚度大于8m范围内，采用650m左右井距，按照油砂体不规则布井，包括老井埕北古4井在内，部署6口定向井作为布井方案一；在此基础上将定向井改为水平井，采用定向井与水平井组合的布井方式，又设计了4种布井方案；再以这5种布井方案为基础，分别改变生产压差、水平井段长度及开发方式，共组合设计了14个方案供数值模拟研究。

3）数值模拟结果分析

（1）生产压差优化

为了分析不同生产压差对开发效果的影响，将方案1～10分成5组，每组的2个方案只是生产压差不同，而布井方式、水平井段长度是相同的。

从数值模拟预测的累积采油量与时间的关系曲线可以看出，每组的2个方案相比，开发初期生产压差3MPa的累积采油量明显高于2MPa的；而到开发后期相差不大，生产压差3.0MPa比2.0MPa平均累积多产油0.25×10⁴t，采出程度提高0.07%。总体来看，生产压差对开发效果影响不大，生产压差3.0MPa略好于2.0MPa。

根据埕岛油田东营组已投产井取得的测压资料统计，平均生产压差3.14MPa。埕北古4井试采期间，6～12mm油嘴，生产压差为2.7～6.2MPa，因此，埕北古4井区定向井生产压差取3.0MPa。

（2）布井方式优化

在生产压差优选的基础上，对5种不同布井方式进行了优化研究。

从数值模拟预测的5种不同布井方案的累积采油量与时间的关系曲线可以看出，定向井与水平井组合的方案比纯定向井方案开发效果好些，累积采油量最大差值为2.8×10⁴t；井多的方案比井少的方案好，但井数越多总投资也越高，因此，需要对各方案进行经济评价^[6]

胜利石油管理局建设项目经济技术评估咨询公司.胜利海上埕岛油田1999年产能建设方案（经济评价报告）.1999.，以进一步优选布井方式。

根据经济评价结果分析，布井方案3的经济效益最好，即埕北古4井区采用3口定向井，1口水平井的布井方式。

（3）水平井段长度优化

优选方案（即方案3：水平井段500m、3口定向井1口水平井、生产压差3MPa）以后，改变水平井段长度，进行模拟计算，研究水平井段对开发效果的影响。根据数值模拟结果分析，随着水平井段长度的增加，单位长度累积采油量增加幅度减小。结合目前钻井工艺水平，水平井段取500m左右较为合适。

开发方式优化 埕北古4井区天然能量比较充足，但地饱压差比相似区块小，开采期间随着地层压力下降油层脱气严重，需注水补充能量开发。又由于本井区距主体已开发区较远，储量规模较小，不宜上大规模的注水设备。具体部署时根据油井钻遇油层及投产情况，考虑构造低部位一口井适时就地取水，补充能量开发。并对枯竭开采和补充能量开采的开发效果进行了对比分析。

根据数值模拟结果分析，注水补充能量的比枯竭开采的累积多产油26×10⁴t，采出程度提高6%。由于注水补充能量的投资较枯竭开采的高，最后推荐枯竭开采和注水补充能量开采两种开采方案进行经济评估。

3.油藏工程方案部署

图1埕北古4井区开发方案井位部署图

根据储集层预测、数值模拟研究、初步经济评价结果，结合相似区块的开发经验，在有井钻遇的落实砂体上部署1口水平井、3口定向井，其中利用埕北古4井老井1口、定向井G4A-1井兼探新砂体Ⅰ、定向井G4A-2井兼探新砂体Ⅱ。新砂体Ⅰ若是满含油则在其上部署1口水平井G4A－平2井。如此埕北古4井区2个落实砂体、2个预测砂体，共部署5口开发井，预计建成年产能力18×10⁴t（图1）。

由于本区只完钻1口井，资料较少，对有井钻遇的主力油砂体油水界面深度以及新砂体是否含油难以准确判断；另外，两种深层储集层预测方法在海上是第一次应用，对其储集层预测精度没有十分把握，由此部署的方案存在较大风险。为了尽最大可能规避风险，要求钻井顺序必须严格按照方案要求实施。

五、结论

利用两种预测方法所描述的储集层，在平面上的分布基本重合，且形状、砂层厚度比较接近。说明这两种方法所描述的储集层基本可信。

埕北古4井区距油田主体远且规模较小，应采用以天然能量开采为主，适时就地取水补充能量的开发模式。

针对埕北古4井区中深层特有的地质和油藏特征，探索和开发了一套储集层预测和开发对策研究的方法和技术，对具有类似地质特征的油田新区具有参考价值。

致谢滩海室周英杰、杜玉山、王军、隋淑玲等高级工程师参加了该项目的部分研究；滩海室范崇海、张强、王海虹、柳文秀、曲全工等同志参加了研究工作。在此一并表示感谢。

主要参考文献

[1]俞寿朋.高分辨率地震勘探.北京：石油工业出版社，1993.

[2]刘企英.利用地震信息进行油气预测.北京：石油工业出版社，1994.

[3]李福垲.黑油和组分模型的应用.北京：科学出版社，1996.

[4]陈钦雷.油田开发设计与分析基础.北京：石油工业出版社，1984.

[5]陈元千.实用油气藏工程方法.东营：石油大学出版社，1993.

[6]中国石油天然气总公司计划局，中国石油天然气总公司规划设计总院编.石油工业建设项目经济评价方法与参数（第二版）.北京：石油工业出版社，1994.

‘捌’ 三维承压水非稳定流模型遗传反演方法

11.3.1 数学模型

含水层参数识别方法

其中，H 为地下水水头[L]；S 为贮水系数[量纲一]；K 为渗透系数[L/T]，二秩对称张量；（x，y，z）为笛卡儿坐标系的坐标；t 为时间；G 为不含边界的研究区域；Γ1 为第一类边界；Γ2 为第二类边界；Γ=Γ1+Γ2 为边界；=G+Γ为包括边界的研究区域；n 为边界单位外法向量；H₀ （x，y，z）为初始水头分布；H₁ （x，y，z，t）为第一类边界上的水头分布；q（x，y，z，t）为第二类边界单位面积上的流量分布。

含水层参数识别方法

其中，N_W抽水井（或注水井）数目；N_L分层数；Q_jk为j号井中k层的出水量。

11.3.2 模型的有限元解法

（1）化边值问题（11-18）为变分问题

若固定 t，并使 H（x，y，z，t）∈C²，则对任意η（x，y，z）∈则有：

含水层参数识别方法

利用奥-高公式，并注意边界条件得：

含水层参数识别方法

（2）剖分求解区域

三维水流问题的有限元法是二维水流问题有限元法的自然推广，所以它的基本思想与解题步骤与二维问题有限元法类似。对于三维问题而言，最简单的单元是四面体，我们现在讨论三维地下水水流问题在剖分为四面体单元时的有限元法。

把区域G剖分为N_e个四面体，四面体的顶点为节点，节点编号i=1，2，3，…，N_P；其中内部节点N₀个，第一边界节点N₁个，第二边界节点N₂个。第i号节点的坐标记为（x_i，y_i，z_i），四面体单元用◆β（β=1，2，…，N_e）表示。

（3）构造基函数和线性形状函数

设单元◆β的四个节点为i，j，k，m。相应的坐标为（x_i，y_i，z_i），…，（x_m，y_m，z_m），四个顶点的水头值H_i，H_j，H_k，H_m。恰好在单元内可确定一线性插值。

含水层参数识别方法

其中：

含水层参数识别方法

其中，V_β为四面体◆β的体积。

L_r可写作如下形式：

含水层参数识别方法

其中系数可通过式（11-19）展开求得。

节点i的基函数定义为：

含水层参数识别方法

我们取形状函数为线性元，则在节点i上水头的分片线性函数可表示为：

含水层参数识别方法

取在边界Γ₁上为0的任意分片光滑函数η为

含水层参数识别方法

式中，m=N₀+N₂，r_i为任意常数。

（4）建立有限元方程

将式（11-22）和式（11-23）代入（11-17）式并考虑r_i的任意性，得有限元方程：

含水层参数识别方法

对初始条件可表示为：

含水层参数识别方法

用全隐式求解常微分方程组（7-24）得隐格式

含水层参数识别方法

式中，=H（x_i，y_i，z_i，t_n），Δt=t_n+1-t_n，T_j=G_j-F_j其中：

含水层参数识别方法

其中，N_t表示以i和j为公共节点的单元个数，K_xxβ，K_yyβ，K_zzβ，S_β为K，S在单元β内的值。

含水层参数识别方法

方程组（11-30）最后可化为AH=B的形式，A为大型稀疏对称矩阵，可用压缩方法存储并求解。如果未知节点太多，用直接解法速度太慢，可用高斯-塞德尔迭代方法求解。

11.3.3 遗传反演方法

三维地下水问题的遗传反演方法和二维地下水问题的反演方法没有本质的差别，所不同的是在二维承压含水层中有导水系数的概念，在三维地下水模型中只有渗透系数的概念，研究区域和外部的水量交换是通过边界条件来实现的，其反演速度更慢，问题可能更复杂。我们通过例子说明三维遗传反演方法。

例 三维承压地下水模型遗传反演问题。承压含水层系统有三个水平层，每层的厚度为20 m。首先作正演计算时从上到下含水层的渗透系数分别为K₁=5 m/d，K₂=2 m/d，K₃=1 m/d，S₁=S₂=S₃=0.0005。从平面上看，含水系统为一正方形，和前述二维承压地下水水流模型算例相同（见图11-2），边长a为1200 m，东、西部边界BC和AD为两条水头均为H₁的河流，切割三个含水层，南、北部边界AB和CD为不透水边界，在含水层的中心有一抽水井，从上层顶端和下层底部抽水。含水系统的上、下均为不透水层。定解问题描述如同方程组（11-18）。

此问题只能用数值方法进行求解，我们采用有限元，并且用直接解法。假设顶部抽水水量Q₁=5000.0 m³/d，底部抽水水量Q₂=2000.0 m³/d。如果A为坐标系的原点，AB为X轴的正方向，AD为Y轴的正反向，含水层系统从底部沿A点向上为Z轴的正方向，则各观测点的资料见表11-7和表11-8。以这些资料为依据反演含水层系统的渗透系数和贮水系数。

表11-7 各观测点水位观测资料一览表（m）

表11-8 各观测点位置一览表（m）

因前面讨论的遗传算法及其组合有8个，本三维模型计算采用直接方法解方程组，速度较慢。我们采用前面比较好的方法即优体克隆+子体优生遗传算法进行反演计算，其计算时的参数分别为：渗透系数的初始取值区间为（0，100.0 m/d），贮水系数S的初始取值区间为（0，0.1），遗传代数Num-Gen=1000，种群数Pop-Size=50，交叉概率P_c=0.7，变异概率P_m=0.3，评价函数中的alph=0.05。计算结果见图11-5和表11-8。从计算结果看，进化到1000代时目标函数并没有达到稳定状态还有下降的趋势。此时各参数已和真值很接近。结果比较见表11-9。反演的结果比较好（表11-10）。

表11-9 优体克隆+子体优生遗传算法反演结果表

表11-10 反演参数和真值比较表

图11-5 优体克隆+子体优生遗传反演方法进化代数与目标函数计算结果图

‘玖’ 两功率分配策略和三功率分配策略的区别

本文主要研究的是多用户MIMO系统的功率分配策略。从单用户系统开始，分析了其“注水算法”，把其与“均匀分配”策略作了比较；进而研究了其推广算法——上行多用户“迭代注水”算法，针对其缺少用户间功率分配导致用户较多时出现的频谱效率降低的现象提出了一种改进的算法，加入了用户间的功率分配内容；进而结合时分策略对改进的算法进行了进一步的简化；并通过建立与下行多用户系统对偶的上行多用户系统的方式，把改进的上行多用户“迭代注水”算法做了相应的改动，运用到下行，来取得下行多用户系统的可达速率，对下行广播信道的容量区域这一开放性问题作了初...