地下储气库（以下简称储气库）运行过程中，由于注采周期交替频繁，储气库压力一直处于动态变化中^[1-4]，稳定的压力对于储气库的安全运行至关重要，而实时监测储气库压力并不现实，尤其是储层压力和井底压力。已有学者围绕储气库压力计算进行了相关研究，如于本福等人^[5]给出了储气库储层压力分布预测方法，唐立根等人^[6]研究了储气库井底流入动态，岳三琪等人^[7]提出了储气库井底压力计算方法。而储气库不同阶段压力之间是相互关联、相关影响的，若仅计算分析储层压力或井底压力，则难以全面掌握储气库压力的动态变化，因此亟需建立一种储气库一体化压力计算方法。

储气库压力主要包括储层压力、井底压力和井口压力，因此本文涉及的压力计算方法，即储层压力、井底压力和井口压力的计算方法。储气库注气过程中，压力传输过程为由井口通过井筒传输至井底，再由井底通过储层传输至远端储层。压力计算流程为：1）利用储层压力计算方法求出储气库的储层压力；2）采用井底压力计算方法，将井底压力与储层压力联系起来，计算得到井底压力；3）利用井口压力计算方法，将井口压力与井底压力联系起来，计算得到井口压力，从而建立储层、井底、井口压力的联系。随着压力不断变化，地层物性参数往往也随之发生改变，如渗透率等^[8-10]；而在注采井注气过程中，受注气量、注气时间等因素的影响，探测半径等也可能在不断变化，如将上述参数视为某一常数，则与储气库的实际生产情况不符。获得这些物性参数的常见方法是进行定期或不定期试井测试^[11]，而试井测试往往需要关井一段时间，影响储气库的正常生产。储气库运行过程中，通常会监测井口压力，因此可以通过综合前述储气库压力的关系，依据现场的静、动态资料，选择一种合适的优化算法来拟合井口压力，使计算出的井口压力与实测井口压力达到最优拟合，从而得到储层压力和井底压力。

常见的优化算法较多，如粒子群优化算法^[12]、遗传算法^[13]及梯度下降法^[14]等。其中，粒子群优化算法因收敛速度快，全局寻优能力强，被广泛应用于求解多目标优化问题和非线性规划问题^[15]。基本粒子群优化算法的核心思想是把种群粒子的位置抽象为目标优化问题的解，通过综合分析并学习粒子和种群最佳位置，动态调整种群中各粒子的位置和速度，经过多次迭代，逼近空间中的最优解^[16]。由于随着迭代次数增加，基本粒子群优化算法的收敛速度逐渐变慢、种群多样性变差，因此，笔者通过引入自动计算惯性因子的策略和变异策略，增加了种群的多样性，使算法的收敛速度加快，得到了一种改进粒子群优化算法。笔者在前人研究的基础上^[17-25]，基于改进粒子群优化算法，结合储层压力、井底压力和井口压力的计算方法，建立了储气库注气过程中一体化压力及地层参数计算方法，通过智能化拟合井口压力计算得到储气库压力，获得合理的地层物性参数，从而指导储气库的安全生产。

1.   压力计算方法

1.1   储层压力

假设储气库注采井储层为圆形封闭空间，井筒注气过程可以看作流体呈平面径向流流入地层，则注气过程的拟压力方程^[26]为：

初始条件为：

内边界条件为：

外边界条件为：

式中：$\psi (·)$为拟压力函数；$p$为t时刻任一位置r处的压力，MPa；r为储层中任意一点与井轴的距离，m；η为导压系数，m²/s；t为注气时间，s；p_i为储层初始压力，MPa；r_w为井筒有效半径，m；G_g为注气重力流量，N/s；$\mu$为黏度，mPa·s；K为渗透率，mD；h为储层有效厚度，m；r_e为封闭储层外边界半径，m；$\phi$为孔隙度；Z 为井筒内气体偏差因子；$C\left( p \right)$为气体等温压缩系数，MPa⁻¹。

在前述边界及初始条件下，将导压系数线性化，取$C(p){\text{ }} = {\text{ }}C({p_{\rm{i}}}){\text{ }} = {\text{ }}1/{p_{\rm{i}}}$，则导压系数简化为η=Kp_i/(ϕμ)；联立式（1）—式（6），对其进行求解得：

结合拟压力函数定义，若将${{\mu {{Z}}} \mathord{\left/ {\vphantom {{\mu {{Z}}} p}} \right. } p}$视为常数，令${{{C_1}{\text{ = }}\mu {{Z}}} \mathord{\left/ {\vphantom {{{C_1}{\text{ = }}\mu {\text{Z}}} p}} \right. } p}$，则有：

将重力流量转换成体积流量，并结合真实气体状态方程，则有：

式中：${q_\text{g} }$为气井注气量，10⁴m³/d；T_f为地层温度，K；p_a为标准大气压，MPa；$\overline {Z}$为地层温度、地层平均压力下的天然气偏差因子；T_a为标准状况温度，K。

将式（8）和式（9）代入式（7），得到注气后储层压力计算公式为：

式中：${p_r}$为t时刻探测半径r处的储层压力，MPa。

1.2   井底压力

通过产能方程建立储层压力与井底压力（井底流压）之间的联系，依据1.1节求得的储层压力，结合注气量，从而计算得到井底压力。产能方程一般分为二项式产能方程和指数式产能方程，二项式产能方程为：

式中：${p_{{\text{wf}}}}$为井底压力，MPa；A为二项式产能方程层流系数，$\text{MPa} /({\text{1}}{{\text{0}}^{\text{4}}}{{\text{m}}^3} \cdot {{\text{d}}^{{{ - 1}}}})$；B为二项式产能方程紊流系数，$\text{MPa} /{({\text{1}}{{\text{0}}^{\text{4}}}{{\text{m}}^3} \cdot {{\text{d}}^{{{ - 1}}}})^2}$。

指数式产能方程为：

式中：C为产能方程系数，${\text{1}}{{\text{0}}^{\text{4}}}\left({{\text{m}}^3} \cdot {{\text{d}}^{{{ - 1}}}} \right) \mathord{\left/ {\vphantom {{\left( {{\text{1}}{{\text{0}}^{\text{4}}}{{\text{m}}^3} \cdot {{\text{d}}^{{{ - 1}}}}} \right)} {\text{MP}{\text{a}^{2n}}}}} \right. } {\text{MP}{\text{a}^{2n}}}$；n为渗流指数。

1.3   井口压力

注采井注气不一定是连续的，可能存在临时关井的情况。因此，井口压力计算方法可分为注气时井口压力计算方法和关井时井口压力计算方法。

1.3.1   注气时井口压力

注采井注气期间，天然气可看作干气，由于注采井井筒长度相对其直径非常大，井筒中的垂直流动可视为以井筒距离为自变量的一维流动，即拟稳定状态流动，满足井筒压力变化规律。根据Cullender和Smith提出的井筒内流体流动的垂直管流方程^[27]，建立井底压力与井口压力之间的关系，结合1.2节求得的井底压力，进而推导得到注气时的井口压力。研究过程中，假设储气库垂直注采井井筒内的流动为单相均质流体的稳定流动，而对于单相均质流体的稳定流动过程，考虑气体在管内与注采井筒壁面所产生的摩擦压力损失，依据注采井筒内的能量平衡，利用伯努利方程，得到以注采井井底为基准面，通过垂直井筒任一流动截面上单位体积气体的压力、动能、位能和由于摩擦产生的压力损失之和为一常数，即：

式中：p为任一流动截面上单位体积气体的压力，MPa；${\rho _{\text{g}}}$为任一流动截面上气体的密度，kg/m³；g为重力加速度，m/s²；d为注采井井筒内径，m；v为井筒内流体的流动速度，m/s；$\Delta {p_\text{f}}$为任一流动截面上单位体积气体的摩擦压力损失，MPa；C₀为井筒内的总压力损失常数，MPa。

式中：f为井筒摩擦系数；L为注采井筒井底深度，m。

将式（14）代入式（13），并微分得：

若忽略流体动能的影响，则以注采井井底为基准面，垂直井筒任一流动截面上压力变化的微分方程为：

其中

式中：Re为雷诺数；$\varepsilon$为注采井井筒内管道的粗糙度，m。

根据气体流动方程，管内流动气体的密度${\rho _{\text{g}}}$为：

式中：γ_g为井筒内气体的相对密度；R为摩尔气体常数，$\text{MPa} \cdot {\text{m}^3} \cdot {\left( {\text{kmol} \cdot \text{K}} \right)^{ - 1}}$；T为井筒内气体温度，K。

井筒内任意压力、温度状态下的气体流速为：

式中：p_sc为拟对比压力，MPa；T_sc为拟对比温度，K。

将式（17）—式（19）代入式（16），整理得：

将式（20）中的T和Z取平均值$\overline{T}、{\overline{Z}}_\text{f}$，并对其从井底到井口进行积分，并整理得到井口压力计算公式为：

式中：p_h为井口压力，MPa；$\overline T$为井筒的平均温度，K；${\overline Z _\text{f}}$为井筒的平均压缩因子。

1.3.2   关井时井口压力

由于关井时气体不流动，则井口压力为：

将式（18）代入式（23），可进行求解。

上述2种井口压力计算方法均采用迭代法和逐点计算法进行求解。

2.   改进的粒子群优化算法

2.1   基本粒子群优化算法

基本粒子群优化算法的主要步骤为：

1）产生初始种群。随机产生N个粒子作为初始种群，各粒子为需要拟合的参数。

式中：${X_{id}}$表示第i个粒子第d维的位置；${X_{i\text{min}}}$表示第i个粒子第d维最小值的位置；${X_{i\text{max}}}$表示第i个粒子第d维最大值的位置；${F_1}$为产生$\left( {0,1} \right)$区间上随机数的函数。

2）计算目标函数值。目标函数值的大小表明粒子的优劣程度，以计算井口压力与实测井口压力的最小均方根误差为目标函数，计算公式为：

式中：${p_{\text{h} j}}$为第j天的计算井口压力，MPa；$p_{\text{h} j}^{\text{real} }$为第j天的实测井口压力，MPa；m为注气时间，d。

3）粒子速度和位置的更新公式。

速度更新公式为：

位置更新公式为：

式中：$\omega$为惯性因子；${v_{id}}$表示第i个粒子第d维的速度；${P_{id}}$表示第i个粒子第d维的位置；${P_{{g} d}}$表示全局最优解第d维的位置。

惯性因子为非负值，其值较大时，全局寻优能力强，局部寻优能力弱；较小时，全局寻优能力弱，局部寻优能力强。

2.2   改进粒子群优化算法

改进粒子群优化算法的主要步骤为：

1）产生初始种群。与基本粒子群优化算法一致，各粒子为需要拟合的参数，包括平均渗透率、探测半径等地层参数，同式（24）。

2）计算目标函数值。将产生的初始种群，作为初始条件代入式（10），计算出储气库的储层压力，根据式（11）或（12）计算出注采井井底压力，将其代入式（21）或（23），计算得到注采井的井口压力，然后再通过式（25）计算得到各粒子的目标函数值。

3）更新粒子的速度和位置。通过引入自动调整$\omega$的策略，可以对全局寻优性能和局部寻优性能进行调整，从而达到加速收敛的效果，节约时间代价。其中，粒子的速度更新公式同式（26），位置更新公式同式（27）。

式中：${f_i}$为第i个粒子目标函数值；${f_{\min }}$为目标函数值的最小值；$\overline f$为目标函数值的平均值。

4）变异策略。引入变异策略，随机选择第i个粒子的第d维进行变异，增加种群的多样性。

式中：${F_2}$表示产生$\left( { - 1,1} \right)$区间上随机数的函数。

重复步骤2）～4），直至种群中粒子表现满足收敛结束条件时停止。

此时种群的全局最优粒子，即为使计算井口压力与实测井口压力达到最优拟合的平均渗透率、探测半径等地层参数，从而求得储层压力和井底压力，进一步基于确定的地层参数预测后期的压力。

基于改进粒子群优化算法的压力和地层参数计算流程如图1所示。

图  1  基于改进粒子群优化算法的压力和地层参数计算流程

Figure  1.  Flow chart of pressure and formation parameter calculation based on improved PSO algorithm

下载: 全尺寸图片 幻灯片

3.   实例验证

呼图壁气藏位于准噶尔盆地南缘山前褶皱带第三排构造带的东端，储层岩性主要为细砂岩、不等粒砂岩和粉砂岩，孔隙以原生粒间孔为主，为受岩性构造控制、带边底水的中孔中渗砂岩贫凝析气藏^[28]。呼图壁储气库由呼图壁气田开发中后期改建而成，于2013年投入注采运行，目前处于第9个注采周期。以呼图壁储气库的3口注采井来验证模型，其中实测井口压力为井口油压。3口注采井的基础参数见表1，其中，X1井渗透率为2017年和2018年3次注气期间通过试井解释获得，最小为8.52 mD，最大为11.26 mD，平均为9.97 mD；X2井渗透率为2018年和2020年3次注气期间通过试井解释获得，最小为9.31 mD，最大为16.93 mD，平均为12.75 mD。采用基本粒子群优化算法和改进粒子群优化算法，计算了注采井X1井不同运行时间下的井口油压和地层参数；采用改进粒子群优化算法，计算了注采井X2井和X3井不同运行时间下的井口油压和地层参数，并分析了计算井口油压与实测井口油压的相关性，结果见表2和图2。

表  1  3口实例井的基础参数

Table  1.  Basic parameters of three example wells

井名	井深/m	井筒有效 半径/m	孔隙度	天然气 相对密度	天然气黏度/ （mPa·s）	产能方程 类型	产能系数A(C)	产能系数B(n)
X1井	3 529.00	0.076	0.165	0.78	0.020	二项式	0.004 3	0.030 4
X2井	3 553.75	0.076	0.155	0.65	0.020	二项式	0.386 2	0.031 9
X3井	3 582.00	0.088	0.209	0.75	0.015	指数式	1.391 4	0.840 8

下载: 导出CSV 

| 显示表格

表  2  粒子群优化算法结果对比及地层参数

Table  2.  Comparison of results obtained from PSO algorithm and formation parameters

井名	算法	平均渗透率/mD		探测 半径/m	相对误差①，%	决定系数②	运行 时间/s
		注气前期	注气后期
X1井	基本粒子群优化算法	10.48	15.01	305.59	0.12	0.988 9	1 657.35
	改进粒子群优化算法	10.48	15.01	305.59	0.12	0.988 9	1 445.51
X2井	改进粒子群优化算法	11.25	19.99	294.28	0.24	0.989 3	857.23
X3井	改进粒子群优化算法	10.00	29.99	301.26	0.11	0.978 4	418.35
注：①和②分别为计算井口油压和实测井口油压的相对误差和决定系数。

下载: 导出CSV 

| 显示表格

图  2  3口注采井实测与计算井口压力的相关性和压力计算及预测结果

Figure  2.  Correlation between calculated and measured wellhead pressure and results of calculated and predicted pressure of three injection and production wells

下载: 全尺寸图片 幻灯片

从表2和图2可以看出：基本粒子群优化算法和改进粒子群优化算法拟合结果一致，说明改进粒子群优化算法具有很好的稳定性；X1井、X2井和X3井的计算井口油压与实测井口油压的相对误差分别为0.12%，0.24%和0.11%，误差相对较小，且决定系数均大于0.97，两者的相关性较强，证明计算井口油压与实测井口油压基本一致，说明拟合效果较好，因而储层压力和井底压力的计算结果可信度较高；计算出的渗透率与试井解释的渗透率基本一致，从而在一定意义上验证了计算出地层参数的合理性，同时也说明随着注气的进行，导致储气库压力不断变化，进而影响地层参数；X1井、X2井和X3井最后10 d的预测井口油压与实测井口油压拟合得也较好，说明预测的井口油压也具备一定的参考价值。

4.   结　论

1）综合储层压力、井底压力和井口压力的计算方法，可以进行地下储气库注气井一体化压力计算。

2）基于改进粒子群优化算法，建立了一种地下储气库注气过程中一体化压力及地层参数计算方法，通过智能拟合计算井口压力与实测井口压力，可以获得地下储气库井底压力和储层压力，进而得到地层参数，并进行压力预测。

3）以呼图壁储气库的3口注采井为例，验证了地下储气库注气过程中一体化压力及地层参数计算方法的可靠性和计算精度。