钻遇窄安全压力窗口地层时，易出现气侵等井下异常工况，严重威胁着钻井安全^[1]。控压钻井（managed pressure drilling，MPD）技术通过精确监测出口流量并控制井筒压力，达到抑制气侵并提升井筒安全的目的^[2]。但井侵气体通常运移至井筒中部甚至是中上部才开始剧烈膨胀，因此早期检测气侵对于MPD系统安全有效抑制气侵至关重要。随着故障诊断技术在各领域的广泛应用，气侵检测在硬件诊断的基础上逐渐融入统计检测、最优估计等理论^[3]。另外，利用数据驱动方法进行气侵诊断的研究工作正逐步兴起^[4]，但其训练速度、可解释性及泛化能力需要进一步探究。同时国内外学者已经开始结合观测器理论来提升气侵早期检测效果。Zhou等^[5]建立包含MPD节流阀的简单水力学模型，设计观测器并通过估计井底流量的变化判断气侵是否发生。Hauge等^[6]通过应用自适应观测器实现了对于地层与井筒间流量进出的定量描述，但其研究过程依赖产能指数已知这一假设。Zhang等^[7]建立参数化描述的井筒模型并设计对应的非线性观测器，通过实时监测调整参数的变化趋势来实现井筒异常工况的早期识别，但该方法中将模型时变项集中到调整参数上可能会给参数反演带来额外的不确定性。徐宝昌等^[8-9]通过坐标变换方法构造自适应观测器，以实现井下气侵工况参数和状态的联合估计，但所建两相流模型中气体膨胀率未能实时更新，且根据特殊坐标变换方法构造自适应观测器的条件十分严苛，难以适用于其他非线性系统。因此笔者建立以“双形”气泡表示法为基础的简化两相流模型，并基于可解耦不确定参数的匹配条件设计自适应观测器，用于实时估计井下气侵工况的不确定参数与未知状态，在此基础上结合霍特林T²检验统计构建气侵检测方法。

1 简化气液两相流模型

实际钻井过程在气侵工况下为气液固三相流动，有学者^[10]通过考虑固相岩屑对环空压降及流动参数的影响来建立井筒三相流模型，其往往具备较高的模型精度，但随之带来了更高的复杂度和解算成本，无法保证方法在现场应用的实时性和鲁棒性。事实上，大多数基于模型的控制和估计方法都限制了它们所应用模型的结构和复杂度^[11]。现有研究常使用漂移流模型建立井筒两相流系统^[12]，用于描述气侵过程伴随的环空复杂流动压力演变规律，然而这些研究更多关注钻井参数的瞬态分析，鲜有学者在其基础上研究控制和估计问题。

为模拟侵入井筒内气体的运动，将整个井筒看作一个固定边界的控制体积（control volume，CV），且CV中仅包含一个气泡，气泡的描述采用“双形”气泡表示法^[13]。如图1所示，气泡的形状以两种不同的方式表示，而选取何种气泡则取决于所要描述井筒气侵工况下的物理特性。实际气侵过程中，井筒内每个气泡在形状和物理性质上都可用泰勒气泡（Taylor bubble，TB）描述，而为描述气体在井筒中的纵向分布，引入拉伸气泡（stretched bubble，SB）。

井筒中液相的一维流动质量守恒可以表示为

式中，m_L为液体质量，kg；ρ_Li和ρ_Lo分别为入口和出口处的液体密度，kg/m³；Q_Li和Q_Lo分别为入口和出口处的液体体积流量，m³/s。

利用质量与密度的关系将导数展开，可得

式中，${\stackrel{-}{\rho }}_{\mathrm{L}}$为平均液体密度，kg/m³；V_L为液体体积，m³，气侵发生阶段与CV的体积近似相等。

假设等温条件，引入液体的线性状态方程：

式中，P_L为环空井底压力，Pa；P₀为参考压力，Pa；ρ₀为标称密度，kg/m³；β为液体体积模量，Pa。

对式（3）求导，可得如下关系：

通过结合式（2）和式（4），可得

式中，$\frac{\mathrm{d}{V}_{\mathrm{L}}}{\mathrm{d}t}$由两种效应产生，即气侵出现时气体占据井筒体积的变化与钻进过程中钻柱所占体积的变化。

因此可表示为

式中，V_G为井侵气体体积，m³；V_D为钻柱在CV中所占体积，m³。

通过将式（6）代入式（5），将气体膨胀项与低阶液体模型耦合，可得

侵入井筒内的气体质量守恒为

式中，m_G为气体质量，kg；q_Gi为从储层侵入井筒内的气体质量流量，kg/s；q_Go为 MPD系统循环排气的质量流量，kg/s。

利用质量与密度关系展开导数并重组方程，可得

式中，ρ_G为气体密度，kg/m³。

钻遇高压气藏时，若井底压力小于地层压力，则地层中的未知气体会侵入井筒并在环空中向上运移。针对此过程，引入线性产能公式^[14]：

式中，J_G为气藏产能指数，kg/（Pa·s）；P_R为地层压力，Pa。

用式（10）替代式（9）中的气侵质量流量，可得两相流模型气体膨胀项：

此外，忽略气液之间的摩擦力，建立两相混合动量平衡方程^[15]：

其中

$\tau \left(v\right)=-\frac{8\mu D_{\mathrm{h}}}{D_{\mathrm{L}}^2}\left(v+K_{\mathrm{C}}{v}_{\mathrm{D}}\right)$

式中，Q_M为井口出口流量，m³/s；m_M为气液混合质量，kg；A_CV为环空横截面积，m²；P_C为井口回压，Pa；L_Ci为环空内外周长和，m；L_CV为井筒长度，m；z_TBb为气泡底部位置，m；z_TBt为气泡顶部位置，m；v_Li为入口处的液体流速，m/s；v_M为气液混合流速，m/s；g为重力加速度，m/s²；θ为井筒与垂直方向夹角，rad；τ（v）为剪切应力，N/m²； D_h为水力直径，m；D_L为兰姆直径，m；μ为钻井液黏度，Pa·s；K_C为黏附系数；v为液体流速，m/s；v_D为钻速，m/s。

选取SB以描述气体实际运移规律，当气侵正在发生且侵入点固定时，SB底部位置不变；当气侵停止时，气泡底部位置变化率为气体运移速度。表示为

其中

$v_{\mathrm{G}}=\gamma v_{\mathrm{s}}+(1-\gamma )v_{\mathrm{b}}.$

式中，z_SBb为SB底部位置，m；v_SBb为SB底部速度，m/s；v_G为气体运移速度^[16]，m/s； v_s为段塞流态气体运移速度，m/s；v_b为泡状流态气体运移速度，m/s；γ为含气率。

气泡顶部的位置则与气体运移和气体膨胀相关：

式中，z_SBt为SB顶部位置，m；α为SB占据环空横截面积的比例。

综上，简化气液两相流模型方程总结如表1所示。

该模型虽然仅由一组显式的常微分方程组成，但仍可不失精度地模拟气侵工况下气液两相的动态流动行为，便于后续在此基础上设计自适应观测器。

2 基于自适应观测器的气侵检测方法

实际气侵过程中气藏产能指数J_G未知，且无法准确判断侵入何种气体，因此可设不确定参数θ₁=J_G/ρ_G，θ₂=（dρ_G/dt）/ρ_G。针对气侵过程，令x₁=P_L，x₂=V_G，x₃=Q_M，简化两相流模型式（7）、（11）、（12）可以改写为状态空间形式：

其中

$\begin{array}{c}{C}_1=\frac{Q_{\mathrm{L}\mathrm{i}}\left(\beta -P_0\right)}{V_{\mathrm{C}\mathrm{V}}},C_2=\frac{Q_{\mathrm{L}\mathrm{i}}}{V_{\mathrm{C}\mathrm{V}}},C_3=\frac{\beta {\rho }_{\mathrm{L}0}}{{\rho }_0{V}_{\mathrm{C}\mathrm{V}}},\\ C_4=\frac{P_{\mathrm{R}}\left(\beta {\rho }_0-P_0{\rho }_0+\beta {\rho }_{\mathrm{L}0}\right)}{2{\rho }_0{V}_{\mathrm{C}\mathrm{V}}},C_5=\frac{P_{\mathrm{R}}{\rho }_0-\beta {\rho }_0+P_0{\rho }_0-\beta {\rho }_{\mathrm{L}0}}{2{\rho }_0{V}_{\mathrm{C}\mathrm{V}}},\\ C_6=\frac{\beta {\rho }_0-P_0{\rho }_0+\beta {\rho }_{\mathrm{L}0}}{2{\rho }_0{V}_{\mathrm{C}\mathrm{V}}},C_7=C_8=\frac{1}{2V_{\mathrm{C}\mathrm{V}}},\\ C_9=P_{\mathrm{R}},C_{10}=\frac{L_{\mathrm{C}\mathrm{i}}{L}_{\mathrm{C}\mathrm{V}}\tau \left(v_{\mathrm{M}}\right)-A_{\mathrm{C}\mathrm{V}}{P}_{\mathrm{C}}-{\rho }_0{V}_{\mathrm{C}\mathrm{V}}g\mathrm{c}\mathrm{o}\mathrm{s}\theta }{{\rho }_0{L}_{\mathrm{C}\mathrm{V}}},\\ C_{11}=\frac{A_{\mathrm{C}\mathrm{V}}}{{\rho }_0{L}_{\mathrm{C}\mathrm{V}}},C_{12}=\frac{\alpha {\rho }_0{A}_{\mathrm{C}\mathrm{V}}g\mathrm{c}\mathrm{o}\mathrm{s}\theta -L_{\mathrm{C}\mathrm{i}}\tau \left(v_{\mathrm{M}}\right)}{\alpha {\rho }_0{V}_{\mathrm{C}\mathrm{V}}}.\end{array}$

式中，C₁~C₁₂均为常数。

2.1 自适应观测器设计

针对式（15）模型，考虑如下形式的非线性系统：

式中，$\mathit{x}\in {\mathbb{R}}^n$为状态向量；$\mathit{y}\in {\mathbb{R}}^m$和$\mathit{y}\in {\mathbb{R}}^p$分别表征输入和输出；$\mathit{\theta }\in {\mathbb{R}}^r$为不确定参数; Ω（y，u）: ${\mathbb{R}}^p\times {\mathbb{R}}^m⟶{\mathbb{R}}^n$为已知向量; Φ（x）: ${\mathbb{R}}^n\times {\mathbb{R}}^m⟶{\mathbb{R}}^s$和Ψ（x）: ${\mathbb{R}}^n\to {\mathbb{R}}^{q\times r}$为未知的非线性函数向量；$\mathit{A}\in {\mathbb{R}}^{n\times n}，\mathit{B}\in {\mathbb{R}}^{n\times s}，\mathit{E}\in {\mathbb{R}}^{n\times q}$和$\mathit{C}\in {\mathbb{R}}^{p\times n}$为已知常数矩阵。

此外，不失一般性地假设rank（C）=p且rank（E）=q。

控制输入u是一致连续且有界的，在保证状态x有界的条件下，对（C，A）是可观测的，且q≤p。此外，对系统（16）提出假设。

假设1   函数Ω（y，u）在其定义域上是连续的，Φ（x）和Ψ（x）关于x满足Lipschitz条件，存在已知正实数γ₁和γ₂使下式成立：

假设2   不确定参数θ有界，即存在已知正实数γ₃使下式成立：

当回归量依赖于未知状态变量时，假设2是设计自适应观测器的标准假设。在应用中，可以通过利用θ的实际性质来设置边界γ₃。

假设3   EΨ（x）满足持续激励（persistent excitation，PE）条件，即存在正常数T、k₁和k₂使得对所有t≥0都有

假设3确保了对于系统不确定参数估计的准确性。然而在实际应用中，由于存在外部干扰等系统未知输入，很难验证回归量是否满足PE条件。

传统观测器匹配条件^[17]为

满足上述匹配条件的实质是能够在系统输出y的一阶导数中完全观测到不确定参数向量。但传统观测器匹配条件在非线性系统的适用性上相对受限，文献^[18]已经验证结构更为简单的MPD单相流水力学系统结构并不能满足传统观测器匹配条件。

针对系统（16），设计自适应状态观测器时基于以下结构假设：

假设4   矩阵A，B，C和E满足

式（21）意味着从输出到Φ（x）和Ψ（x）θ的相对阶至少为2，式（22）为改进的观测器匹配条件。为进一步解释式（21）、（22）的含义，从相对阶的角度来分析，首先对输出y求一阶导数并将式（21）代入可得

再对输出y求二阶导数可得

由式（23）、（24）可知，假设4表明不确定项不直接存在于可测输出y中，需要对输出y求两次导数才能得到关于Φ（x）和Ψ（x）θ的显式表达式，这意味着能够从更微小的输出变化中观测到完整的系统结构。

当假设4成立时，自适应观测器的设计基于以下引理：

引理1   rank（CAE）=rank（E）当且仅当存在矩阵P=P^T＞0，H以及G满足如下关系时成立：

当系统（16）满足假设1~4时，考虑如下形式的类Luenberger观测器^[18]：

不确定参数估计向量$\widehat{\mathit{\theta }}$是以下自适应律的解：

式中，W为不确定参数θ估计值的中间变量; $\mathit{L}\in {\mathbb{R}}^{n\times p}$为观测器增益; H与G为待定的常数矩阵; Γ=Γ^T＞0为学习率矩阵。

自适应观测器式（26）~（28）渐近稳定的条件在下面的定理中给出。

定理1   考虑满足假设1~4的系统（16），如果存在正实数ε₁、ε₂和矩阵P=P^T＞0、H、G使得式（25）及下式成立：

$(\mathit{A}-\mathit{L}\mathit{C}{)}^{\mathrm{T}}\mathit{P}+\mathit{P}(\mathit{A}-\mathit{L}\mathit{C})+{\epsilon }_1\mathit{P}\mathit{B}{\mathit{B}}^{\mathrm{T}}\mathit{P}+{\epsilon }_2\mathit{P}\mathit{E}{\mathit{E}}^{\mathrm{T}}\mathit{P}+$

那么不确定参数估计值$\widehat{\mathit{\theta }}$收敛于参数真实值θ，状态估计值$\widehat{\mathit{x}}$渐近收敛至真实值x。在文献^[18]中，对上述定理进行了详细的推导证明，并通过Lyapunov稳定性理论证明了估计误差的渐近稳定性。式（25）为观测器匹配条件下的矩阵等式，保证了在输出的二阶导数中能够观测到全部的不确定参数。式（29）表示Lyapunov函数的导数小于0，保证了观测器的收敛性。

通过在等式约束式（25）下求解式（29）来计算观测器增益相对困难，解决这类问题较为常见的方法是将其转化为线性矩阵不等式（linear matrix inequalities，LMI）优化问题。首先，令L=P^-1M，其中$\mathit{M}\in {\mathbb{R}}^{n\times p}$。其次，分别利用Schur补性质^[19]和文献^[20]中的理论将观测器增益的求解问题转化为如下形式的LMI优化问题：

其中

$\mathit{\Lambda }={\mathit{A}}^{\mathrm{T}}\mathit{P}+\mathit{P}\mathit{A}-{\mathit{C}}^{\mathrm{T}}{\mathit{M}}^{\mathrm{T}}-\mathit{M}\mathit{C}+{\epsilon }_1^{-1}{\gamma }_1^2{\mathit{I}}_n+{\epsilon }_2^{-1}{\gamma }_2^2{\gamma }_3^2{\mathit{I}}_n.$

式中，η为一个正的标量。当η越趋近于0时，线性矩阵不等式（31）越接近于观测器匹配条件下的矩阵等式（25）。因此需要在上述优化问题的求解中使η最小化，以求取最优的观测器增益H、G和L。

2.2 气侵检测方法

采用霍特林T²检验统计^[21]完善气侵检测方法，该方法由数据驱动离线设计与在线气侵检测两部分组成。

假设有两组数据：未发生气侵时的离线控压数据y_i（i=1，···，N）和在线测量控压数据y_k。根据离线控压过程数据计算均值与方差的无偏估计：

定义检验统计量为

若一个可接受的误报率为α，则相应的阈值为

式中，F_α（m，N-m）是自由度为m和N-m的F分布。

在线气侵检测过程中，根据以下逻辑判断气侵是否发生：

基于自适应观测器的气侵检测方法可以总结如图2所示。

3 仿真分析

3.1 模型验证

OLGA全动态多相流模拟器作为当前世界领先的非稳态多相流模拟软件^[22]，计算结果被世界各大石油公司所认可。将建立的简化两相流模型与OLGA仿真数据对比，以验证该模型应用在气侵早期检测中的有效性。参考现场某井的垂直井段，设置模型参数：L_CV=4000 m，L_Ci=1.237 m，A_CV=0.043 m²，θ=0 rad，P₀=1.01×10⁵ Pa，ρ₀=1198.3 kg/m³，β=1.66×10⁹ Pa，μ=0.05 Pa·s，ρ_G=311.5 kg/m³，P_R=5×10⁷ Pa，J_G=5×10^-6 kg/（Pa·s），P_C=1.2×10⁶ Pa，Q_Li=0.01 m³/s，D_h=0.14 m，D_L=0.05 m，v_D=0.001 m/s，K_C=0.05，α=0.9。设置初始条件：P_L（0）=4.9×10⁷ Pa，V_G（0）=0 m³，Q_M（0）=0.01 m³/s，z_SBb（0）=z_SBt（0）=0 m。

图3为简化两相流模型与OLGA仿真的环空井底压力和井口出口流量对比。其中，前400 s为未发生气侵的单相流动过程；400~600 s为气侵过程（红色区域）；600 s后气侵停止，气泡开始向上运移。

仿真结果表明：在气侵过程中，简化两相流模型计算结果与OLGA仿真数据匹配良好，模型计算精度能够满足气侵早期检测的需要。

3.2 气侵检测方法仿真实验

针对简化两相流模型式（15），应用本文自适应观测器进行仿真分析，并采用OLGA模拟器生成数据作为仿真实验的真实值数据来源。初始状态设置为P_L（0）=4.9×107 Pa，${\widehat{P}}_{\mathrm{L}}（0）$=4.85×10⁷ Pa，${\widehat{V}}_{\mathrm{G}}（0）$=0 m³，Q_M（0）=0.01 m³/s，${\widehat{Q}}_{\mathrm{M}}（0）$=0.013 m³/s；不确定参数初值设置为${\widehat{\theta }}_1（0）$=2.3×10^-8，${\widehat{\theta }}_2（0）$=-3.75×10^-4。

考虑假设1、2，设计Lipschitz常数为γ₁=0，γ₂=97.5，γ₃=4×10^-4。同时，结合实际数据及经验分析，选取学习率矩阵为，再取ε₁=ε₂=50。最后，利用式（30）、（31）计算观测器增益矩阵得

$\begin{array}{c}L=\left[\begin{array}{cc}0.29& -4.81\times {10}^6\\ -8.71\times {10}^{-8}& 3.54\\ -1.18\times {10}^{-8}& 1.326\end{array}\right],\\ G=\left[\begin{array}{cc}-8.27\times {10}^{-14}& 0\\ 1.35\times {10}^{-8}& 0\end{array}\right],\\ H=\left[\begin{array}{cc}-1.78\times {10}^{-13}& 1.73\times {10}^{-6}\\ 8.83\times {10}^{-7}& -1.5\end{array}\right].\end{array}$

仿真结果如图4、5所示。图4（a）为两个不确定参数的估计值和真实值，图4（b）为两个不确定参数的估计误差，图5（a）为井底压力、气体体积、出口流量3个状态的观测器估计值和真实值（OLGA数据），图5（b）为3个状态的估计误差。

由图4可知，不确定参数均能在60 s内收敛至真实值，收敛速度快且估计结果准确。图5仿真结果表明：在气侵工况下，本文自适应观测器能够使估计值在60 s内跟踪上正在变化的井底压力和出口流量，且后续的估计误差几乎为0，体现了观测器良好的跟踪性能。

由式（15）可知，不确定参数对气体体积的估计有着极大影响，自适应观测器在60 s时对气体体积的估计误差为0.21 m³，在120 s时估计误差缩小至0.15 m³，在180 s后估计误差几乎为0。表明本文自适应观测器在不确定参数给定初值与真实值相差较大时，仍能将不确定参数收敛至真实值。自适应观测器对井筒内气体体积的估计效果良好，估计误差能够保持在很小的区间内并持续缩小至0。

Stokka^[23]、Speers等^[24]和Maus等^[25]已经分别验证了声波法、出入口流量差法和总池体积法在气侵早期检测中的有效性，因此选用以上3种方法进行对比实验。为更好地验证本文方法应用于实际工况的可行性，考虑当钻井液密度存在3%的误差时进行仿真实验，同时将塔里木油田钻井现场采集到的入口流量测量噪声（方差为10^-5的白噪声）加入到仿真计算中。在400 s时地层压力突变为5×10⁷ Pa，用于模拟气侵发生，仿真结果如图6、7所示（红色区域表示气侵过程）。

如图6所示，自适应观测器能够在气侵发生后迅速跟踪上井底压力和出口流量的变化，并能够实时、准确地估计出侵入井筒内的气体体积，且观测过程几乎不受噪声扰动影响。图7表示气侵工况下声波传播时间、出入口流量差和总池体积的变化。

为合理验证本文气侵检测方法的优势并结合现场实际工况分析，4种方法的可接受误报率α均设置为0.02。根据式（35），检验阈值分别计算得到$J_{\mathrm{t}\mathrm{h}，T^2}$（本文气侵检测方法）=10.21，$J_{\mathrm{t}\mathrm{h}，T^2}$（声波法）= $J_{\mathrm{t}\mathrm{h}，T^2}$（出入口流量差法）= $J_{\mathrm{t}\mathrm{h}，T^2}$（总池体积法）=5.53。针对4种方法分别计算检验统计量，对比结果如图8所示（红色区域表示气侵过程）。

图8仿真结果表明：本文气侵检测方法基于自适应观测器对井下状态的准确估计，能够使检验统计量快速超过阈值，在5 s内检测出气侵。对比其他3种方法的气侵检测时间，该方法早于声波法（检测效率提高71.4%）和出入口流量差法（检测效率提高80.0%），且相较于总池体积法表现出显著的性能提升（检测效率提高92.3%）。对比其他3种方法的检验统计量变化显著程度，该方法远大于其他检测方法，可有效避免误报警。

此外，为充分说明本文气侵检测方法的有效性，对比分析了不同气侵程度下4种检测方法的气侵检测时间。由式（10）可知，气侵程度主要由钻遇的地层压力决定。因此在仿真分析中，假定当前井底压力为4.9×10⁷ Pa，设置钻遇的地层压力分别为4.92×10⁷、4.94×10⁷、4.96×10⁷、4.98×10⁷ 和5×10⁷ Pa，以表示5种不同大小的气侵程度。如图9所示，4种方法的气侵检测时间都随气侵程度的减小而增大，但声波法、出入口流量差法和总池体积法在气侵程度较小时的检测效果更差，在钻遇4.92×10⁷ Pa的地层压力时，无法在气侵发生过程中及时预警。对比分析可知，特别在气侵程度较小的情况下，本文气侵检测方法相较于声波法、出入口流量差法和总池体积法表现出更好的诊断性能，能够及时且灵敏地检测出气侵。

4 结论

（1）基于“双形”气泡表示法建立的简化气液两相流模型，通过与OLGA仿真数据对比，验证了该模型在复杂度适中的情况下能保证良好的计算精度。

（2）设计的基于观测器匹配条件的自适应观测器，用于实时估计井下气侵工况的不确定参数（气藏产能指数及气体密度）与未知状态（气体体积），结合霍特林T2检验统计构建了气侵检测方法。

（3）当气侵正在发生时，新的自适应观测器能够实时估计侵入井筒的气体体积，并准确跟踪井底压力和出口流量的变化，气侵检测方法相较于声波法、出入口流量差法和总池体积法，可以在不同气侵程度下快速且显著地检测出气侵，能够为抑制气侵赢得时间。

（4）井底压力在气侵发生时会呈现“先增后减”的变化趋势，导致本文中气侵检测方法的检验统计量会在后期出现阶段性下降；在未来的研究中可对检验统计方法进行合理地改进以解决这一问题。

参考文献