随着油气资源勘探不断深入，深层和超深层油气资源已成为我国油气勘探开发的重要领域^[1–5]。然而，油气环境日益复杂，导致深井、超深井钻探时面临着诸多技术挑战，如地层条件复杂、高温高压、施工周期长等，对安全和高效地钻进提出了严格要求^[6–10]。在这一背景下，压井技术作为保障井下作业安全和高效进行的关键手段备受关注，而传统的压井方法主要依赖于固定的立压和套压变化曲线进行控制，缺乏对井筒内复杂多相流型、流体相态和流动规律的实时响应能力。为了解决这一问题，范军等人^[11–12]考虑了储层渗透率和钻开储层厚度等因素的影响，建立了动态多相流动方程组。孙宝江等人^[13–14]针对深水井控的特点，考虑深水钻井外部的多温度梯度环境和天然气水合物相变，建立了七组分井筒多相流控制方程。O. L. A. Santos^[15]建立的深水井涌模型考虑了气液间的滑脱及空隙率的影响，但未考虑井眼形状变化。H. V. Nickens^[16]建立的多相流模型考虑气液相质量变化及混合物的动量变化，并考虑了两相滑脱、井眼形状等对压力的影响。J. O. L. Nunes等人^[17]建立的解析模型主要适用于计算稳态时的情况，主要采用段塞流流型进行计算，在气侵量较大时较为适用。虽然已有许多学者将影响井筒气液混合流体物性参数及运移规律的关键因素纳入考虑范围，但由于地质条件和井筒流体流动特性的复杂性，现有模型仍然存在较大的计算误差，无法准确描述井筒内的气液分布状态和压力变化规律。

目前，随着传感技术、物联网技术的发展和智能化数字钻井技术的进步，石油钻采监测数据量剧增，且呈现工业大数据特点^[18–23]。其中，井下随钻数据、井口实测数据等多源信息不仅反映了复杂环境下井筒流体物性参数和运动参数的实时变化，还包括所钻地层岩石性质及压力体系分布等信息。因此，笔者针对传统模型存在的计算误差大的问题，结合井口实测数据及井下温压数据等多源信息，研究不确定性流动参数变化与实时数据间的映射关系，提出了实时数据约束下不确定性参数的智能反演方法；针对压井初期少量数据约束下模型的多解性问题，构建了基于连续时间序列的模型全局训练优化方法与动态随机种群训练优化方法，实现了对井筒内气液分布状态的精确表征、压力变化规律的实时预测和多源数据约束下压井曲线的实时更新，提高了压井的安全性和效率。

1.   智能压井模型及其多解性分析

1.1   井筒不确定性参数计算方法

分析多相流动模型计算过程发现，误差主要来自漂移流动模型中气相分布系数C₀与气体滑脱速度v_gm计算方程和多相流摩阻系数f_tp计算方程的不准确性。其中，C₀反映了环空截面处气体分布特征，影响着环空多相流流型转化，并决定了溢流及循环压井过程中环空压力分布的变化幅度；v_gm反映了侵入气体的运移速度分布，决定了气体前缘位置及环空压力分布的变化速度；f_tp反映了裸眼井段及套管内壁对环空流体的阻尼特性，从而影响了环空压力分布的变化幅度。因此，C₀与v_gm的精确性可保证井筒内气液分布状态的精确计算，f_tp的精确性可保证压力变化规律的精确计算。

现有计算C₀与v_gm的方程，一方面依赖于室内试验^[24–27]，无法完全复刻井下复杂温压环境及侵入流体的高压物性，使得漂移流动模型无法针对特定的井筒环境精确计算出空隙率分布；另一方面，由于储层区块复杂的地层特性分布，侵入流体类型具有不确定性（例如，侵入流体为酸性气体时，需考虑其在钻井液中的溶解特性；深层钻井大温差环境条件下，需考虑侵入流体的相变特性等），导致固化的漂移流动模型难以准确预测井筒空隙率分布的变化规律^[27–28]。因此，为了实现井筒空隙率分布和压力变化的精确求解，建立了计算漂移流动关键参数C₀，v_gm及摩阻系数f_tp的修正方程：

式中：C₀为气相分布系数；$Cr_{C_0}$为气相分布修正系数；v_gm为气体滑脱速度，m/s；$Cr_{v_{\mathrm{gm}}}$为气体滑脱速度修正系数；f_tp为两相摩阻系数；$Cr_{f_{\mathrm{tp}}}$为多相流摩阻修正系数；${f_{{C_0}}}$（°），${f_{{{{v}}_{{\rm{gm}}}}}}$（°）和${f_{{f_{{\rm{tp}}}}}}$（°）分别为包含所需计算参数的C₀，v_gm和f_tp的原有经验方程，不同学者建立的经验方程有所差异；ρ_g，ρ_l分别为气相、液相密度，kg/m³；μ_l为液相黏度，Pa·s；${\bar v_{\text{g}}}$为气相平均速度，m/s；D为管道水力直径，m；v_sg，v_sl分别为气相、液相表观速度，m/s；m_g，m_l分别为气相、液相质量流量，kg/s；θ为管道倾角，（°）；E_g为截面含气率；ε为管道水力粗糙度，m。

式（1）—式（3）表明，C₀，v_gm及f_tp的修正方程是原有经验方程分别乘以修正系数$Cr_{C_0}$，$Cr_{v_{\mathrm{gm}}}$及$Cr_{f_{\mathrm{tp}}}$，当经验方程计算误差较大时，可通过调整修正系数保证计算的精确性。因此，引入修正方程，可以将多相流动模型的校正问题转化为多源实时数据约束下多元修正系数的反演求解工作。结合井下随钻数据、井口实测数据等信息，可设定目标函数F：

式中：N为监测点总数；p_m,i(t)为监测点i位置处于t时刻的井筒压力测量值，Pa；p(i, t)为对应监测点i位置处的井筒压力模拟值，Pa。

多源信息约束下多相流动模型的实时修正可具体描述为寻找合适的$Cr_{C_0}$，$Cr_{v_{\mathrm{gm}}}$和$Cr_{f_{\mathrm{tp}}}$，使得目标函数取得极小值，其本质为目标函数F约束下的多目标优化问题。研究采用遗传算法计算F极小值点对应的$Cr_{C_0}$，$Cr_{v_{\mathrm{gm}}}$和$Cr_{f_{\mathrm{tp}}}$，其将随机的参数取值模拟为染色体种群，基于种群适应度的选择和重组过程产生下一代的继承群体，即更适应目标函数的$Cr_{C_0}$，$Cr_{v_{\mathrm{gm}}}$和$Cr_{f_{\mathrm{tp}}}$。重组过程中，选择亲本染色体并重组其遗传物质，产生的子染色体随后流入继承群体，算法的迭代过程即为一系列的连续世代进化过程，其中染色体的平均适应度不断提高，直至达到算法终止条件。其流程主要包括染色体编码，目标函数及染色体选择、重组和进化方案^{[27, 29]}。详细智能压井模型及修正参数优化过程可参考文献[27]。

1.2   智能压井模型多解性分析

由上述分析可知，$Cr_{v_{\mathrm{gm}}}$主要控制气相的运移速度，对井筒压力分布的影响较小；而$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$分别通过影响井筒静液柱压力和环空摩阻实现井筒压力分布的调控，因此在同一目标函数F（$Cr_{C_0}$，Cr_vgm，$Cr_{f_{\mathrm{tp}}}$）约束下，可能存在不同的$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$组合满足F的全局最小化，即$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$对井筒压力的双重调节作用可能导致目标函数F的多解性。

以我国西部某案例井为例^[27,30]进行计算分析。该井井深3 954 m，地温梯度3.0 ℃/100m，钻井排量0.03 m³/s；用水基钻井液，密度1 200 kg/m³，比热容1 800 J/(kg·℃)，导热系数0.82 W/(m·℃)；打开储层厚度10 m，储层比热容832 J/(kg·℃)；钻杆/套管导热系数50 W/(m·℃)；水泥环比热容900 J/(kg·℃)，导热系数0.75 W/(m·℃)。该井采用四级井身结构，三开采用ϕ311.1 mm钻头钻至3 370 m处，下入壁厚11.0 mm的ϕ244.5 mm技术套管，四开采用ϕ215.9 mm钻头配合ϕ127.0 mm钻杆钻至3 954 m。钻具组合为：ϕ215.9 mm钻头+ϕ152.4 mm钻铤×1根+ϕ209.5 mm稳定器+ϕ152.4 mm钻铤×1根+ϕ209.5 mm稳定器+ϕ152.4 mm钻铤×5根+ϕ127.0 mm加重钻杆×12根+ϕ127.0 mm钻杆。

结合该案例井数据，得到不同$Cr_{f_{\mathrm{tp}}}$取值条件下$Cr_{C_0}$变化对井底压力计算结果的影响规律，见图1。从图1可以看出，$Cr_{C_0}$变化使得井底压力呈非线性变化趋势，随着$Cr_{C_0}$增大，井底压力的增长幅度逐渐减小；随着$Cr_{f_{\mathrm{tp}}}$增大，井底压力呈线性增大。从图1中红色虚线可以得出，不同的$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$组合可以产生相同的井底压力。

图  1  不同修正系数组合对井底压力的调控规律

Figure  1.  Regulation law of different correction coefficient combinations on bottom hole pressure

下载: 全尺寸图片 幻灯片

修正系数取值过大或过小，都会导致多相流模型求解出现不收敛的现象。经过对案例井的反复测试求解，设置修正系数$Cr_{C_0}$取值范围为0.50～1.43，$Cr_{f_{\mathrm{tp}}}$取值范围为0.55～1.50，以井底压力实测值与模拟值误差的绝对值为目标函数，绘制不同$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$组合进化前后的误差分布云图，见图2。从图2可以看出，随机初始种群经过遗传算法计算得到的最优解收敛至图2（b）中所示的红色包络线区域，即$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$组合形成的二维空间存在井底压力的“等高线”，表明在目标函数F的约束下，由于少量数据条件下多相流动模型无法准确反演井筒空隙率分布及气相运移规律，从而存在无穷多组$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$组合满足其局部最优化。因此，应综合考虑压井过程中井筒气液流动参数的变化规律，利用实时信息构成的连续时间序列进行多相流动模型的训练与完善。

图  2  单一目标函数约束下修正系数组合的多解性

Figure  2.  Multi-solution of correction coefficient combinations under single objective function constraint

下载: 全尺寸图片 幻灯片

2.   模型全局训练优化方法

2.1   历史节点空隙率对参数计算的影响

通过分析井筒多相流模型求解过程中的空隙率计算流程及其内在原理，可以解决目标函数F的多解性问题。井筒单位空间节点（i，j）空隙率计算流程为：假设空隙率E_g(i，j)，结合连续性方程计算气相真实速度v_g1，结合漂移流动模型计算气相真实速度v_g2；循环迭代E_g，使v_g1与v_g2间的相对误差满足许可值，从而求解空隙率。气相、液相连续性方程为：

式中：A为环空横截面积，m²；v_g，v_l分别为气相、液相真实速度，m/s；t为时间，s；s为空间距离，m；E_l为持液率。

采用一阶向后差分对时间偏导，并采用一阶迎风差分对空间偏导，可以得到当气相、液相速度均为正时j时刻气相、液相连续性方程的离散格式和漂移流动关系式：

式中：t_step为时间步长，s；x_step为空间步长，m。

综合式（7）—式（9），可得气相速度的表达式：

式中：v_g1，v_g2分别为连续性方程、漂移流动模型计算得到的气相真实速度，m/s。

修正系数$Cr_{C_0}$和$Cr_{v_{\mathrm{gm}}}$是通过调控漂移流动模型所计算的气相真实速度v_g2，以实现对空隙率的调控。式（10）所示的连续性方程表明，历史时间节点空隙率E_g(i，j−1)影响气相真实速度v_g1的计算，因此同样对当前时间节点的空隙率计算有显著的影响；结合式（11）所示的漂移流动模型可得，历史时间节点空隙率E_g(i，j−1)同样影响当前时间节点漂移流动模型对气相真实速度v_g2的计算。

历史时间节点空隙率E_g(i，j-1)变化对当前节点连续性方程和漂移模型计算结果的影响规律如图3所示。从图3（a）可以看出，在不同E_g(i，j−1)取值条件下，基于连续性方程计算的气相真实速度v_g1变化幅度较大，而基于漂移模型计算的气相真实速度v_g2变化不明显。同时，E_g(i，j−1)增大，导致当前节点基于连续性方程的气相真实速度的计算曲线随空隙率的变化速率整体增大，曲线形态整体变“陡峭”，表明此时相同的$Cr_{C_0}$变化将对应更大幅度的空隙率波动（见图3（b）），即历史节点空隙率增大，增强了$Cr_{C_0}$的调控能力。而$Cr_{C_0}$对空隙率的大幅调控，又导致下一个节点处空隙率计算结果的变化，因此两者间的耦合促进作用将导致最优解空间形态件下，基于连续性方程计算的气相真实速度v_g1变化幅度较大，而基于漂移模型计算的气相真实速度峭”，表明此时相同的$Cr_{C_0}$变化将对应更大幅度的空隙率波动（见图3（b）），即历史节点空隙率增大，的变化。压井过程中井筒空隙率是随着气体段塞的运移而逐渐增大的，这表明可以通过联立不同时间节点处目标函数F关于$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$的分布规律，消除修正系数组合的多解性。

图  3  历史节点空隙率变化对当前节点连续方程和漂移模型计算结果的影响

Figure  3.  Influence of void fraction change in historical nodes on calculation results of continuity equation of current nodes and drift model

下载: 全尺寸图片 幻灯片

2.2   解空间的形态演变规律

以1.2节所示案例井为例，模拟压井过程中不同时间节点解空间的形态演变规律，结果如图4所示。分析图4中包络线内最优解区域形态的演变过程可知，随着压井过程中井筒空隙率增大，$Cr_{C_0}$的变化幅度逐渐减弱，表明$Cr_{C_0}$的小幅波动导致的井筒静液柱压力变化即可抵消$Cr_{f_{\mathrm{tp}}}$变化导致的环空摩阻变化，因此出现包络线形态逐渐“平缓”的趋势（对应图4（a）—图4（d）中的深蓝色区域）。

图  4  压井过程中解空间形态的演化规律

Figure  4.  Evolution law of solution spatial morphology during well killing

下载: 全尺寸图片 幻灯片

2.3   全局训练下模型最优解收敛过程

综合上述分析可以得出，在实时信息时间序列约束下最优解的搜寻过程实质上是多相流动模型关于井筒气液分布规律的深度学习过程，解空间形态的变化与收敛域的收缩则对应着模型的训练与优化。

结合2.2节分析可得，$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$的最优解组合需满足图4（a）—图4（d）中对应的目标函数群的全局最小化，只需重合图4（a）—图4（d）中包络线区域，对应的收敛区域即为最优解区域。然而，这种图解法难以编程运算，因此将图4（a）—图4（d）中的误差分布叠加至图4（a），叠加结果如图5所示。从图5可以看出，叠加后的红色包络线区域（即计算结果误差值等高线包裹区域）逐渐缩小，对应图5（d）即为全局最优解。

图  5  全局训练优化方法下修正系数的收敛趋势

Figure  5.  Convergence trend graph of correction coefficients under global training optimization method

下载: 全尺寸图片 幻灯片

3.   动态随机种群训练优化方法

全局训练优化方法的引入解决了修正系数组合的多解性问题，且保持了较高的计算精度。然而，修正系数的全局搜索增加了计算耗时，因此需要针对计算条件相对紧张的情况，提出多解性问题的解决方法。

上述压井过程中，计算耗时的增加主要来源于修正系数$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$在每一个计算单元内均保持固定的取值范围，这种静态的分布区间设计忽略了相邻时间节点对应的修正系数组合之间的关联关系，导致了单位计算网格内搜索区间的扩大。因此，基于相邻计算网格对应搜索区间的函数关系，得到了修正系数种群的动态随机种群训练优化方法，该方法以上一时间节点处修正系数收敛值为中心，并选择合适的搜索半径作为修正系数的分布范围，以给定的小范围内井底压力实测值与模拟值误差的局部最小值对应的$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$作为单位网格内计算的修正系数组合，这样的动态搜索区域可以大幅减少单位计算网格内的计算耗时。

动态随机种群训练优化方法下压井初期修正系数的收敛过程如图6所示。从图6可以看出，在给定的小搜索范围区间的约束下，$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$收敛至局部最优解，并以此修正系数组合进行下一时间节点的井筒压力计算，得到对应节点实时数据后，修正系数组合又继续收敛至下一局部最优解，并在实时信息序列的约束下逐步收敛至全局最优解。对比全局训练优化方法解空间形态变化对应的模型训练与优化过程，动态随机种群训练优化方法中修正系数的收敛过程则对应着模型关于井筒空隙率分布及压力变化的学习过程，因此，修正系数由局部最优向全局最优收敛过程的本质是连续实时数据序列约束下模型的训练与优化。

图  6  动态随机种群训练优化方法下修正系数的收敛过程

Figure  6.  Convergence process of correction coefficients for dynamic random population training optimization method

下载: 全尺寸图片 幻灯片

4.   两类方法对比分析及适用性评价

以1.2节所示案例井为例，设定压井时井底压力恒为48.85 MPa，得到全局训练优化方法与动态随机种群训练优化方法下智能压井压力的变化特征。

全局训练优化方法下井筒压力的调控过程如图7所示。从图7可以看出，这种搜索方法可以在极短时间内调控井底压力至设计值，实现了井筒压力的精确调控。与图5对比可以看出，虽然压井初期（压井时间t=150 s之前）的修正系数可以实现井底压力模拟值与实测值的误差最小化，但由于小数据量约束下仍存在修正系数组合的多解性问题，模型仍需进一步训练完善。

图  7  全局训练优化方法对应的智能压井压力变化特征

Figure  7.  Pressure change characteristics of intelligent well killing corresponding to global training optimization method

下载: 全尺寸图片 幻灯片

全局训练优化方法的引入虽然解决了多解性难题，然而实时数据下多相流关键参数的反演计算需要耗费大量的计算资源。以本算例为例，在保证修正系数精度为小数点后一位时，遗传算法所需的初始种群规模为10×10×10，这意味着单位计算网格内多相流动机理与实时数据耦合驱动模型的计算耗时将是单纯依赖模型计算的1 000倍，因此该方法更适用于计算资源充足的情况。

动态随机种群训练优化方法下智能压井曲线的变化特征如图8所示。从图8（b）可以看出，147 s以前，依据实时数据反馈计算得出的井底压力仍小于地层压力。结合图6可以看出，井底压力模拟值与实际值之间的差异来自于动态随机种群训练优化方法导致的搜索范围的缩减，使得遗传算法关于$Cr_{C_0}$，$Cr_{v_{\mathrm{gm}}}$和$Cr_{f_{\mathrm{tp}}}$的多元反演求解收敛至局部最小值。这表明此时多相流模型仍处于训练过程中，无法准确反映井筒空隙率分布，且井底压力也未直接增长至预期值，而呈现出快速增长的趋势。随着采集数据量的增加（147 s后），动态随机种群逐渐转移到全局最优解，$Cr_{C_0}$，$Cr_{v_{\mathrm{gm}}}$和$Cr_{f_{\mathrm{tp}}}$逐渐收敛，表明此时模型对井筒空隙率的学习程度已显著提高，井底压力随时间的延长也逐渐接近地层压力。综上可知，采集的实时信息间接反映了井筒裸眼段岩石物性参数和井筒流体的流动特性，修正后的多相流模型计算精度明显提高，证明了该方法的适用性。

图  8  动态随机种群训练优化方法对应的智能压井曲线变化特征

Figure  8.  Change characteristics of intelligent well killing curve corresponding to dynamic random population training optimization method

下载: 全尺寸图片 幻灯片

综合上述分析可得，修正系数的动态随机种群训练优化方法在解决多解性问题的同时，引发了模型训练期间修正系数组合的收敛过程，使得压井早期井底压力未直接增长至预期值，而是保持迅速增大趋势。与整个压井过程耗时相比，模型训练耗时占比约5.9%，训练过程所需数据量为24组；而在保证相同的种群密度条件下，该方法将单位网格内的计算耗时由10×10×10降低至（r/m）×（r/m）×（r/m）（r为搜索半径，m为计算精度，在本算例中计算耗时为4×4×4），在保证计算精度的前提下极大地减少了计算耗时，拓宽了模型的适用性。

5.   结论与建议

1）多相流动机理与实时数据耦合驱动模型的多解性主要源于$Cr_{C_0}$和$Cr_{f_{\mathrm{tp}}}$对井筒压力的双重调节作用，导致在同一目标函数约束下，可能存在不同的参数组合满足目标函数的局部最优化。

2）历史时间节点空隙率增大，会导致当前节点基于连续性方程的气相真实速度的计算曲线随空隙率的变化速率整体增大，使相同的$Cr_{C_0}$变化对应更大幅度的空隙率波动，从而影响最优解空间形态。

3）分析解空间形态的演变规律发现，可以利用不同节点收敛域重合与累加的方式逐步缩小收敛域范围，实现修正系数组合最优解的全局训练优化。

4）考虑相邻时间节点修正系数组合收敛域的函数关系，发现可以通过缩小遗传算法种群搜索区间的方式，并结合历史时间节点处实时信息序列的联合约束，实现从局部最优解逐步收敛至全局最优解的动态随机种群训练优化。

5）全局训练优化方法可以在极短时间内实现井筒压力的精确调控，但未克服在压井初期小数据量约束下修正系数组合的多解性，同时导致单位网格内计算耗时大幅增长，更适用于计算资源充足的情况。动态随机种群训练优化方法在压井初期以局部最优解代替全局最优解，使得压井早期压力未直接增长至预期值，而是保持迅速增大趋势；但单位网格的计算耗时较小，适用于计算资源紧张的情况。