本发明涉及钻井技术，具体地，涉及一种基于双测点测量的井下气侵工况确定方法及系统。
背景技术：
目前，油气资源勘探开发难度越来越大，很多钻井作业是在以高温高压、窄泥浆密度窗口、多套复杂压力体系等为典型特征的深部地层中进行，导致漏、喷、塌、卡等井下复杂事故频发，钻井效率与安全问题非常突出。而控压钻井技术通过闭环自动控制系统来精确控制整个井眼环空压力分布，使井筒压力控制在安全密度窗口内，可有效节约钻井时间、降低钻井成本、提高钻井作业安全性，目前已成为解决复杂地层窄安全密度窗口问题，保障安全高效钻井的最有效的技术手段之一。
然而由于复杂地层存在着大量的模糊性、随机性和不确定性问题，因此在利用控压钻井作业时，复杂工况是不可能完全避免的。而对井下复杂工况认识不清会致使技术人员决策失误和后期盲目的控制措施，可能导致情况复杂化甚至造成更为严重的事故，其中气侵工况危害最为严重，一旦控制不及时不得当，就可能导致井涌、井喷恶性事故，在此时精确的井下工况解释就显得尤为重要。井下工况解释主要是指依据地面录井数据和井下随钻压力测量数据，结合井筒水力学模型进行实时分析和解释。
目前的解释方法主要是利用井下pwd单点测量数据进行解释，然而依据单点参数测量值具有局限性和不确定性，无法准确反映工况变化。其次，由于井筒流动状态的高度复杂性，导致单纯依靠水力学模型来精确描述井筒流体流动是难以实现的，计算结果与真实值往往存在较大误差。当控压钻井过程中出现气侵复杂工况时，已有的研究技术无法准确的定位气侵层位以及确定气侵量大小等井下参数，直接影响到后续的井筒压力控制过程，导致压力控制滞后与精度低。此外，气侵工况解释理论的缺乏在一定程度上也制约了控压钻井技术的发展。
技术实现要素：
本发明实施例的主要目的在于提供一种基于双测点测量的井下气侵工况确定方法及系统，以实时准确预测井下气侵层位与气侵量。
为了实现上述目的，本发明实施例提供一种基于双测点测量的井下气侵工况确定方法，包括：
步骤1：输入最新历史时刻的状态向量的均值及协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值；其中，最新历史时刻为距目标时刻最近的已知状态向量的均值和已知协方差矩阵的时刻；
步骤2：输入多个状态向量估计值至已创建的井筒水力学正演模型中，得到下一时刻的多个测量向量估计值；
步骤3：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值；
步骤4：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵；
步骤5：判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复步骤1至步骤4，直至下一时刻为目标时刻；
步骤6：根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
在其中一种实施例中，输入最新历史时刻的状态向量的均值及协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值，包括：
根据最新历史时刻的状态向量的均值和最新历史时刻的协方差矩阵，计算最新历史时刻的多个确定采样点的状态向量；
输入多个确定采样点的状态向量至非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值。
在其中一种实施例中，根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值，包括：
根据多个测量向量估计值和下一时刻的测量噪声的方差，计算下一时刻的自协方差测量矩阵；
根据多个状态向量估计值和多个测量向量估计值，计算下一时刻的互协方差矩阵；
根据状态向量估计值、测量向量估计值、自协方差测量矩阵、互协方差矩阵和下一时刻的测量参数，计算下一时刻的状态向量的均值。
在其中一种实施例中，根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵，包括：
根据多个状态向量估计值和最新历史时刻的过程噪声的方差，计算下一时刻的自协方差状态矩阵；
根据自协方差状态矩阵和自协方差测量矩阵，计算下一时刻的协方差矩阵。
在其中一种实施例中，根据多个测量向量估计值和下一时刻的测量噪声的方差，计算下一时刻的自协方差测量矩阵，包括：
根据多个测量向量估计值和多个均值加权权值，计算下一时刻的测量向量估计值的均值；
根据多个测量向量估计值、测量向量估计值的均值、多个协方差加权权值和下一时刻的测量噪声的方差，计算下一时刻的自协方差测量矩阵。
在其中一种实施例中，根据多个状态向量估计值和最新历史时刻的过程噪声的方差，计算下一时刻的自协方差状态矩阵，包括：
根据多个状态向量估计值和多个确定采样点对应的多个均值加权权值，计算下一时刻的状态向量估计值的均值；
根据多个状态向量估计值、状态向量估计值的均值、多个确定采样点对应的多个协方差加权权值和最新历史时刻的过程噪声的方差，计算下一时刻的自协方差状态矩阵。
在其中一种实施例中，根据多个状态向量估计值和多个测量向量估计值，计算下一时刻的互协方差矩阵，包括：
根据多个状态向量估计值、状态向量估计值的均值、多个测量向量估计值和测量向量估计值的均值，计算下一时刻的互协方差矩阵。
在其中一种实施例中，还包括：
根据初始时刻的多个状态向量计算初始时刻的状态向量的均值；
根据初始时刻的多个状态向量和初始时刻的状态向量的均值计算初始时刻的协方差矩阵。
在其中一种实施例中，还包括：
根据下一时刻的第一测点压力值、下一时刻的第二测点压力值和下一时刻的井口的返出流量，确定下一时刻的测量参数。
本发明实施例还提供一种基于双测点测量的井下气侵工况确定系统，包括：
状态向量估计值单元，用于输入最新历史时刻的状态向量的均值及最新历史时刻的协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值；其中，最新历史时刻为距目标时刻最近的已知状态向量的均值和已知协方差矩阵的时刻；
测量向量估计值单元，用于输入多个状态向量估计值至已创建的井筒水力学正演模型中，得到下一时刻的多个测量向量估计值；
状态向量的均值单元，用于根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值；
协方差矩阵单元，用于根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵；
判断单元，用于判断下一时刻是否为目标时刻；
气侵工况确定单元，用于根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
在其中一种实施例中，还包括：
初始单元，用于根据初始时刻的多个状态向量计算初始时刻的状态向量的均值；
根据初始时刻的多个状态向量和初始时刻的状态向量的均值计算初始时刻的协方差矩阵。
在其中一种实施例中，还包括：
测量参数确定模块，用于根据下一时刻的第一测点压力值、下一时刻的第二测点压力值和下一时刻的井口的返出流量，确定下一时刻的测量参数。
本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行计算机程序时实现以下步骤：
步骤1：输入最新历史时刻的状态向量的均值及协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值；其中，最新历史时刻为距目标时刻最近的已知状态向量的均值和已知协方差矩阵的时刻；
步骤2：输入多个状态向量估计值至已创建的井筒水力学正演模型中，得到下一时刻的多个测量向量估计值；
步骤3：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值；
步骤4：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵；
步骤5：判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复步骤1至步骤4，直至下一时刻为目标时刻；
步骤6：根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
本发明实施例还提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现以下步骤：
步骤1：输入最新历史时刻的状态向量的均值及协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值；其中，最新历史时刻为距目标时刻最近的已知状态向量的均值和已知协方差矩阵的时刻；
步骤2：输入多个状态向量估计值至已创建的井筒水力学正演模型中，得到下一时刻的多个测量向量估计值；
步骤3：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值；
步骤4：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵；
步骤5：判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复步骤1至步骤4，直至下一时刻为目标时刻；
步骤6：根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
本发明实施例的基于双测点测量的井下气侵工况确定方法及系统，可以对井下气侵工况进行实时校正计算，实时判别井下气侵复杂工况，准确预测井下气侵层位与气侵量，扩大了井筒水力学的应用范围，提高了计算精度，为后续的井筒压力控制提供重要参考。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明实施例中基于双测点测量的井下气侵工况确定方法的流程图；
图2是本发明实施例中基于双测点测量的井下气侵工况确定方法其中一种实施例的流程图；
图3是本发明实施例中基于双测点测量的井下气侵工况确定方法另一种实施例的流程图；
图4是本发明实施例中基于双测点测量的井下双测点测量装置的示意图；
图5是本发明实施例中基于双测点测量的井下气侵工况确定系统的结构框图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
鉴于目前已有的研究技术无法准确的定位气侵层位以及确定气侵量大小等井下参数，直接影响到后续的井筒压力控制过程，导致压力控制滞后与精度低，本发明实施例提供一种基于双测点测量的井下气侵工况确定方法及系统，可以对井下气侵工况进行实时校正计算，实时判别井下气侵复杂工况，准确预测井下气侵层位与气侵量，扩大了井筒水力学的应用范围，提高了计算精度，为后续的井筒压力控制提供重要参考。以下结合附图对本发明进行详细说明。
图1是本发明实施例中基于双测点测量的井下气侵工况确定方法的流程图。如图1所示，基于双测点测量的井下气侵工况确定方法包括：
s101：输入最新历史时刻的状态向量的均值及协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值；其中，最新历史时刻为距目标时刻最近的已知状态向量的均值和已知协方差矩阵的时刻。
s102：输入多个状态向量估计值至已创建的井筒水力学正演模型中，得到下一时刻的多个测量向量估计值。
s103：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值。
s104：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵。
s105：判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复s101至s104，直至下一时刻为目标时刻；
s106：根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
由图1所示的流程可知，本发明先输入最新历史时刻的状态向量的均值及协方差矩阵至非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值，再输入多个状态向量估计值至井筒水力学正演模型得到下一时刻的多个测量向量估计值，进一步计算得到下一时刻的状态向量的均值和下一时刻的协方差矩阵。然后判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复s101至s104，直至下一时刻为目标时刻，最后根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。利用该基于双测点测量的井下气侵工况确定方法，可以对井下气侵工况进行实时校正计算，实时判别井下气侵复杂工况。
在执行s101之前，还包括：创建非线性系统状态函数模型和井筒水力学正演模型。
一实施例中，s101具体包括：根据最新历史时刻的状态向量的均值和最新历史时刻的协方差矩阵，计算最新历史时刻的多个确定采样点的状态向量。以及，输入多个确定采样点的状态向量至非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值。
具体实施时，可以通过如下公式计算最新历史时刻的多个确定采样点的状态向量：
其中，k-1时刻为最新历史时刻，xi,k-1为k-1时刻的第i个确定采样点的状态向量，
λ＝α2(n+κ)-n；
其中，α为扩展因子，影响确定采样点(sigma点)围绕状态向量的均值的分布状态，一般取值[0,1]；κ为可调参数，对于本发明中的高斯分布，当状态向量为n维时，κ＝3-n。例如，若确定采样点的状态向量为二维状态向量，则κ取值为1。
以及，非线性系统状态函数模型如下：
xi,k/k-1＝f(xi,k-1)i＝0,1,…,2n；
其中，k时刻为最新历史时刻的下一时刻，xi,k/k-1为最新历史时刻的下一时刻的第i个状态向量估计值。
s102中，井筒水力学正演模型如下：
zi,k/k-1＝h(xi,k/k-1)i＝0,1,…,2n；
其中，zi,k/k-1为最新历史时刻的下一时刻的第i个测量向量估计值。
图2是本发明实施例中基于双测点测量的井下气侵工况确定方法其中一种实施例的流程图。如图2所示，实施例中，s103具体包括：
s201：根据多个测量向量估计值和下一时刻的测量噪声的方差，计算下一时刻的自协方差测量矩阵。
s202：根据多个状态向量估计值和多个测量向量估计值，计算下一时刻的互协方差矩阵。
s203：根据状态向量估计值、测量向量估计值、自协方差测量矩阵、互协方差矩阵和下一时刻的测量参数，计算下一时刻的状态向量的均值。
一实施例中，s201包括：根据多个测量向量估计值和多个均值加权权值，计算下一时刻的测量向量估计值的均值。以及，根据多个测量向量估计值、测量向量估计值的均值、多个协方差加权权值和下一时刻的测量噪声的方差，计算下一时刻的自协方差测量矩阵。
具体实施时，可以通过如下公式计算下一时刻的测量向量估计值的均值：
其中，
以及，通过如下公式计算计算下一时刻的自协方差测量矩阵：
其中，
一实施例中，s202包括：根据多个状态向量估计值、状态向量估计值的均值、多个测量向量估计值和测量向量估计值的均值，计算下一时刻的互协方差矩阵，具体公式如下：
其中，
图3是本发明实施例中基于双测点测量的井下气侵工况确定方法另一种实施例的流程图。如图3所示，s104具体包括：
s301：根据多个状态向量估计值和最新历史时刻的过程噪声的方差，计算下一时刻的自协方差状态矩阵。
s302：根据自协方差状态矩阵和自协方差测量矩阵，计算下一时刻的协方差矩阵。
一实施例中，s301包括：根据多个状态向量估计值和多个确定采样点对应的多个均值加权权值，计算下一时刻的状态向量估计值的均值。以及，根据多个状态向量估计值、状态向量估计值的均值、多个确定采样点对应的多个协方差加权权值和最新历史时刻的过程噪声的方差，计算下一时刻的自协方差状态矩阵。
具体实施时，通过如下公式计算下一时刻的状态向量估计值的均值：
其中，
通过如下公式计算下一时刻的自协方差状态矩阵：
其中，pk/k-1为k时刻的自协方差状态矩阵，qk-1为k-1时刻的过程噪声的方差，满足零均值高斯白噪声分布wk-1～n{0,qk-1}，wk-1为k-1时刻的系统过程噪声。
wi(c)为第i个确定采样点对应的协方差加权权值，如下：
其中，α为扩展因子，β为状态向量先验分布信息参数。对于本发明而言，β取高斯分布的最优值2。
一实施例中，本发明还包括：根据k时刻(下一时刻)的互协方差矩阵和k时刻的自协方差测量矩阵，计算k时刻的卡尔曼增益矩阵，如下：
其中，kk为k时刻的卡尔曼增益矩阵。
将上述公式带入s103中，则s103具体包括：根据下一时刻的状态向量估计值的均值、下一时刻的测量向量估计值的均值、下一时刻的卡尔曼增益矩阵和下一时刻的测量参数，计算下一时刻的状态向量的均值，如下：
其中，
将上述公式带入s104中，则s104具体包括：根据下一时刻的自协方差状态矩阵、下一时刻的自协方差测量矩阵和下一时刻的卡尔曼增益矩阵，计算下一时刻的协方差矩阵，如下：
其中，pk为k时刻的协方差矩阵。
一实施例中，本发明还包括：根据初始时刻的多个状态向量计算初始时刻的状态向量的均值；根据初始时刻的多个状态向量和初始时刻的状态向量的均值计算初始时刻的协方差矩阵。具体公式如下：
其中，
一实施例中，本发明还包括：根据下一时刻的第一测点压力值、下一时刻的第二测点压力值和下一时刻的井口的返出流量，确定下一时刻的测量参数。具体公式如下：
zk＝[pm,1(tk),pm,2(tk),qm,s(tk)]t；
其中，zk为下一时刻的测量参数，pm,1(tk)为下一时刻的第一测点压力值，pm,2(tk)为下一时刻的第二测点压力值，qm,s(tk)为下一时刻的井口的返出流量。
图4是本发明实施例中井下双测点测量装置的示意图。如图4所示，井下双测点测量装置包括：泥浆脉冲发生器、两个测量短节、流量计，压力表、节流阀和泥浆泵。泥浆泵用于抽取泥浆池中的泥浆送入井口，压力表用于测量井口的压力，节流阀用于对井口压力进行加压或泄压。两个测量短节安装在井底近钻头处(图4中的测点1和测点2)。测量短节之间采用有线钻杆的方式连接，保证了双测点测量数据的同步传输，最终测量数据通过上方的泥浆脉冲发生器发送至地面，本领域技术人员可通过地面计算机解码获取测点1和测点2的测量参数。考虑到有线钻杆的成本问题，测点1和测点2之间的距离不易过大，一般取值为10m-100m比较合适。需要说明的是，在本发明中测点1和测点2的测量参数均为压力，通过井下双测点测量装置的测点1和测点2可分别获得下一时刻的第一测点压力值和下一时刻的第二测点压力值；通过井下双测点测量装置的流量计可获得下一时刻的井口的返出流量。除压力之外，测点1和测点2还可以测量温度、三轴振动、钻压、扭矩和介电常数等参数。
在实施s106时，此时的下一时刻为目标时刻，且下一时刻的状态向量的均值相当于下一时刻的状态向量xk，即
xk＝[qkick(tk),hkick(tk)]t；
此时，qkick(tk)为目标时刻(下一时刻)的气侵量，hkick(tk)为目标时刻的气侵位置。因此，可以通过下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
在本发明的其中一个具体实施例中，设最新历史时刻为第3时刻，则最新历史时刻的下一时刻为第4时刻；设目标时刻为第6时刻。本发明通过以下步骤确定目标时刻的气侵量和目标时刻的气侵位置。
1、根据第3时刻的状态向量的均值和第3时刻的协方差矩阵，计算多个第3时刻的确定采样点的状态向量：
其中，xi,3为第3时刻的第i个确定采样点的状态向量，
2、输入多个第3时刻的确定采样点的状态向量至非线性系统状态函数模型中，得到多个第3时刻的状态向量估计值：
xi,4/3＝f(xi,3)i＝0,1,…,2n；
其中，xi,4/3为第3时刻的第i个状态向量估计值。
3、根据多个第4时刻的状态向量估计值和多个均值加权权值，计算第4时刻的状态向量估计值的均值：
其中，
4、根据多个第4时刻的状态向量估计值、第4时刻的状态向量估计值的均值、多个协方差加权权值和第3时刻的过程噪声的方差，计算第4时刻的自协方差状态矩阵：
其中，p4/3为第4时刻的自协方差状态矩阵，wi(c)为第i个确定采样点对应的协方差加权权值，q3为第3时刻的过程噪声的方差。
5、输入多个第4时刻的状态向量估计值至井筒水力学正演模型中，得到多个第4时刻的测量向量估计值：
zi,4/3＝h(xi,4/3)i＝0,1,…,2n；
其中，zi,4/3为第4时刻的第i个测量向量估计值。
6、根据多个第4时刻的测量向量估计值和多个均值加权权值，计算第4时刻的测量向量估计值的均值：
其中，
7、根据多个第4时刻的测量向量估计值、第4时刻的测量向量估计值的均值、多个协方差加权权值和第4时刻的测量噪声的方差，计算第4时刻的自协方差测量矩阵：
其中，
8、根据多个第4时刻的状态向量估计值、第4时刻的状态向量估计值的均值、多个第4时刻的测量向量估计值和第4时刻的测量向量估计值的均值，计算第4时刻的互协方差矩阵：
其中，
9、根据第4时刻的互协方差矩阵和第4时刻的自协方差测量矩阵，计算第4时刻的卡尔曼增益矩阵：
其中，k4为第4时刻的卡尔曼增益矩阵。
10、根据第4时刻的状态向量估计值的均值、第4时刻的测量向量估计值的均值、第4时刻的卡尔曼增益矩阵和第4时刻的测量参数，计算第4时刻的状态向量的均值：
其中，
11、根据第4时刻的自协方差状态矩阵、第4时刻的自协方差测量矩阵和第4时刻的卡尔曼增益矩阵，计算第4时刻的协方差矩阵：
其中，p4为第4时刻的协方差矩阵。
12、判断下一时刻是否为目标时刻：第4时刻不是目标时刻，则重复上述1-11，此时最新历史时刻为第4时刻，其下一时刻为第5时刻；第5时刻不是目标时刻，继续重复上述1-11，最新历史时刻变更为第5时刻，其下一时刻为第6时刻；第6时刻是目标时刻，则进行下一步：
13、根据第6时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置：
通过1-13，可以得到目标时刻的气侵量qkick(t6)和目标时刻的气侵位置hkick(t6)。
综上，本发明实施例的基于双测点测量的井下气侵工况确定方法先输入最新历史时刻的状态向量的均值及协方差矩阵至非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值，再输入多个状态向量估计值至井筒水力学正演模型得到下一时刻的多个测量向量估计值，进一步计算得到下一时刻的状态向量的均值和下一时刻的协方差矩阵。然后判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复上述内容，直至下一时刻为目标时刻，最后根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置，可以对井下气侵工况进行实时校正计算，实时判别井下气侵复杂工况。
基于同一发明构思，本发明实施例还提供了一种基于双测点测量的井下气侵工况确定系统，由于该系统解决问题的原理与基于双测点测量的井下气侵工况确定方法相似，因此该系统的实施可以参见方法的实施，重复之处不再赘述。
图5是本发明实施例中基于双测点测量的井下气侵工况确定系统的结构框图。如图5所示，基于双测点测量的井下气侵工况确定系统包括：
状态向量估计值单元，用于输入最新历史时刻的状态向量的均值及最新历史时刻的协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值；其中，最新历史时刻为距目标时刻最近的已知状态向量的均值和已知协方差矩阵的时刻。
测量向量估计值单元，用于输入多个状态向量估计值至已创建的井筒水力学正演模型中，得到下一时刻的多个测量向量估计值。
状态向量的均值单元，用于根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值。
协方差矩阵单元，用于根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵。
判断单元，用于判断下一时刻是否为目标时刻。
气侵工况确定单元，用于根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
在其中一种实施例中，还包括：
初始单元，用于根据初始时刻的多个状态向量计算初始时刻的状态向量的均值；
根据初始时刻的多个状态向量和初始时刻的状态向量的均值计算初始时刻的协方差矩阵。
在其中一种实施例中，还包括：
测量参数确定模块，用于根据下一时刻的第一测点压力值、下一时刻的第二测点压力值和下一时刻的井口的返出流量，确定下一时刻的测量参数。
本发明实施例的基于双测点测量的井下气侵工况确定系统先输入最新历史时刻的状态向量的均值及协方差矩阵至非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值，再输入多个状态向量估计值至井筒水力学正演模型得到下一时刻的多个测量向量估计值，进一步计算得到下一时刻的状态向量的均值和下一时刻的协方差矩阵。然后判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复上述内容，直至下一时刻为目标时刻，最后根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置，可以对井下气侵工况进行实时校正计算，实时判别井下气侵复杂工况。
基于同一发明构思，本发明实施例还提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行计算机程序时实现以下步骤：
步骤1：输入最新历史时刻的状态向量的均值及协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值；其中，最新历史时刻为距目标时刻最近的已知状态向量的均值和已知协方差矩阵的时刻；
步骤2：输入多个状态向量估计值至已创建的井筒水力学正演模型中，得到下一时刻的多个测量向量估计值；
步骤3：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值；
步骤4：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵；
步骤5：判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复上述步骤1至步骤4，直至下一时刻为目标时刻；
步骤6：根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
本发明实施例的计算机设备先输入最新历史时刻的状态向量的均值及协方差矩阵至非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值，再输入多个状态向量估计值至井筒水力学正演模型得到下一时刻的多个测量向量估计值，进一步计算得到下一时刻的状态向量的均值和下一时刻的协方差矩阵。然后判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复上述内容，直至下一时刻为目标时刻，最后根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置，可以对井下气侵工况进行实时校正计算，实时判别井下气侵复杂工况。
基于同一发明构思，本发明实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现以下步骤：
步骤1：输入最新历史时刻的状态向量的均值及协方差矩阵至已创建的非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值；其中，最新历史时刻为距目标时刻最近的已知状态向量的均值和已知协方差矩阵的时刻；
步骤2：输入多个状态向量估计值至已创建的井筒水力学正演模型中，得到下一时刻的多个测量向量估计值；
步骤3：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差、下一时刻的测量噪声的方差和下一时刻的测量参数，计算下一时刻的状态向量的均值；
步骤4：根据多个状态向量估计值、多个测量向量估计值、最新历史时刻的过程噪声的方差和下一时刻的测量噪声的方差，计算下一时刻的协方差矩阵；
步骤5：判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复上述步骤1至步骤4，直至下一时刻为目标时刻；
步骤6：根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置。
本发明实施例的计算机可读存储介质先输入最新历史时刻的状态向量的均值及协方差矩阵至非线性系统状态函数模型中，得到最新历史时刻的下一时刻的多个状态向量估计值，再输入多个状态向量估计值至井筒水力学正演模型得到下一时刻的多个测量向量估计值，进一步计算得到下一时刻的状态向量的均值和下一时刻的协方差矩阵。然后判断下一时刻是否为目标时刻；如果下一时刻不是目标时刻，则重复上述内容，直至下一时刻为目标时刻，最后根据下一时刻的状态向量的均值确定目标时刻的气侵量和目标时刻的气侵位置，可以对井下气侵工况进行实时校正计算，实时判别井下气侵复杂工况。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。