0 引言

油相渗流微分方程是描述油相在地下渗流过程中全部物理现象的统一微分方程或微分方程组, 油相在多孔介质渗流的过程中, 遵循质量守恒定理, 将油相的状态方程、岩石状态方程及其他辅助方程代入到质量守恒方程中则可得到油相渗流的基本微分方程。文献[1]采用Laplace变换求解油相渗流微分方程, 为地下流体渗流微分方程求解引入新的方法。20世纪60年代初, 文献[2]建立双重介质达西渗流的数学模型, 文献[3]对双重介质达西渗流数学模型采用Laplace变换进行求解。20世纪末, 文献[4]建立三重介质油藏数学模型, 并采用Laplace变换进行了求解。文献[5]考虑井储效应和表皮效应的渗流数学模型进行研究, 求解得到Laplace空间解。文献[6]也对这类三重介质做了理论研究。总结前人研究成果，外边界封闭单纯介质和多重介质储层流体渗流过程中，压力波波及到边界之后，采用Laplace变换求解，获得Laplace空间解析解和实空间解析解，都与复合Bessel方程的零点有关。Van Everdinger和W.Hurst研究成果中只给出了部分参数复合Bessel方程前3个零点，在前人的研究基础之上，结合现代优化算法对复合Bessel方程零点求解方法及分布规律进行探索。

1 基本模型及解析

油相在地下渗流基本模型的假设条件:流体在储层中进行单相渗流; 储层中流动的油相是微可压缩; 储层为均质水平等厚; 油相在储层中的流动服从达西定律; 不考虑储层污染和井储效应。

基本渗流数学模型参见文献[6]。为了简化研究问题, 通常需要把某些物理量进行参数的无量纲化, 为此引入无量纲参数, 对基本的渗流数学模型、边界条件和初始条件进行无量纲化, 然后对时间变量进行Laplace变换, 将偏微分方程转换为二阶常微分方程, 该方程含正则奇点, 在正则奇点邻域进行级数展开并整理可得方程的通解, 结合边界条件和初始条件可得到渗流微分方程的特解^[7]。

有界封闭地层井底定产模型解

$ p\left( {r,t} \right) = {p_{\rm i}} - \frac{{q\mu }}{{2\pi kh}}\left\{ \frac{{2\eta t}}{{(r^2_{\rm De} - 1)r^2_{\rm w}}} + \frac{{r^2_{\rm De}}}{{r^2_{\rm De} - 1}}\left(\ln \frac{{{r_{\rm De}}}}{{{r_D}}} + \frac{{r^2_D}}{{2r^2_{\rm De}}} \right)- \frac{{3r^4_{\rm De} - 4r^2_{\rm De}\ln {r_{\rm De}} - 2r^2_{\rm De} - 1}}{{4{{(r^2_{\rm De} - 1)}^2}}} - \right.\\ \left. \quad \pi \sum\limits_{n = 1}^\infty {\frac{{\exp( - \alpha ^2_n\frac{{\eta t}}{{r^2_{\rm w}}}){J_1}({r_{\rm De}}{\alpha _n})[{Y_1}({\alpha _n}){J_0}({r_D}{\alpha _n}) - {J_1}({\alpha _n}){Y_0}({r_D}{\alpha _n})]}}{{{\alpha _n}[J^2_1({\alpha _n}{r_{\rm De}}) - J^2_1({\alpha _n})]}}} \right\}, $

(1)

式中: p为压力，0.01 MPa；r为半径，cm；t为时间，s；p_i为原始地层压力，0.1 MPa；q为井下产量，cm³/s；μ为黏度，mPa·s；k为渗透率，μm²；h为储层有效厚度，cm；η为导压系数，η=k/(ϕμc_t)，其中ϕ为孔隙度，c_t为综合压缩系数，(0.1 MPa)^-1；r_w为井筒半径，cm；r_De为外边界无量纲半径，即外边界半径与井筒半径之比，r_De=r_e/r_w，其中r_e为油藏外边界半径，cm；α_n为复合Bessel方程的第n个零点；J₀、J₁为第一类零阶、一阶Bessel函数；Y₀、Y₁为第二类零阶、一阶Bessel函数。

有界封闭地层井底定压模型解

$ q\left( t \right) = \frac{{4\pi hk({p_{\rm i}} - {p_{\rm wf}}({r_{\rm w}},t))}}{u}\sum\limits_{n = 1}^\infty {\frac{{\exp( - \frac{{\alpha ^2_nkt}}{{\phi \mu {c_t}r^2_w}})[{Y_1}({r_{\rm De}}{\alpha _n}){J_1}({\alpha _n}) - {Y_1}({\alpha _n}){J_1}({\alpha _n}{r_{\rm De}})]}}{{{r_{\rm De}}[{Y_0}({r_{\rm De}}{\alpha _n}){J_0}({\alpha _n}) - {J_0}({r_{\rm De}}{\alpha _n}){Y_0}({\alpha _n})] + {Y_1}({\alpha _n}){J_1}({r_{\rm De}}{\alpha _n}) - {Y_1}({r_{\rm De}}{\alpha _n}){J_1}({\alpha _n})}}} , $

(2)

其中: p_wf为井底流压, 0.1 MPa。

2 复合Bessel函数的特性及零点求解策略

通过基本渗流微分方程的求解, 得到的复合Bessel方程是第一类零阶、一阶Bessel函数, 第二类零阶、一阶Bessel函数的复合表达形式, 并且与无量纲半径有关。为更好地研究复合Bessel方程Y₁(r_Deα_i)J₀(α_i)-Y₀(α_i)J₁(r_Deα_i)=0, 定义复合Bessel函数

$ f(x,{r_{\rm De}}) = {Y_1}({r_{\rm De}}x){J_0}\left( x \right) - {Y_0}\left( x \right){J_1}({r_{\rm De}}x), $

(3)

其中:x为函数的自变量, x∈(0, +∞]。

通过式(3), 可绘制带参数的复合Bessel函数图形, 见图 1, 复合Bessel函数为多峰函数, 其零点有无穷多个。

在具体问题研究中, 需要用到复合Bessel函数前面有限个有序零点, 结合复合Bessel函数的特征, 为了更好地求解函数的零点, 将函数零点的求解转化为多峰函数的极值求解问题, 可将式(3)变换为

$ F(x,{r_{\rm De}}) = |f(x,{r_{\rm De}})|。$

(4)

3 复合Bessel函数零点求解算法

针对有界地层模型, 具体工程计算中只需要取前有限个有序零点即可, 如何获得这个复杂数学模型的高精度零点一直是工程应用中精确计算的难解之题。因此, 寻找一种好的方法获得复合Bessel函数有限个有序零点值供设计人员选择不同的设计方案尤为重要^[8-10]。为获得多峰函数特定范围内所有峰值解以及提高所有峰值解的精度、提高整个求解效率等问题, 许多学者采用各种智能算法、进化算法以及各种改进的智能算法以提高多峰函数优化解的性能指标, 并取得了许多研究成果^[11-17]。文献[18]为了尽可能寻找多峰函数的全部极值点及提高寻优精度, 在普通粒子群(particle swarm optimization, PSO)和量子粒子群(quantum-behaved particle swarm optimization, QPSO)算法的基础之上, 针对粒子吸引点限定在全局最优位置和个体历史最优位置的超维矩形内, 提出带交叉算子量子粒子群(crossover-operator quantum particle swarm optimization, CQPSO)算法; 文献[19]针对QPSO算法在求解高维多峰函数难以优化粒子每一维和易陷入局部极值点问题, 采取多次蒙特卡洛随机观测值, 进行前后代粒子逐维对比优化, 并构造一种新的调控收缩扩张系数; 文献[20]针对基本遗传算法在求解多峰函数时很难找到全部最优解的不足, 提出一种并行免疫遗传算法求解该问题, 算法引入郭涛算法的多父体杂交思想, 并借鉴小生境机制将算法的优化过程进行分解。上述研究成果所给优化方法全部是解决高维多峰函数全局极值问题, 对于复合Bessel函数这样一维多峰函数仍然使用, 通过图 1可知, 本次求解目标函数F(x, r_De)极小值全部为0, 不存在全局最小值, 在具体工程应用中, 必须求出该函数的前面有限个有序极值点才行, 如何在求解过程中保证不遗漏极值点且精度高才是求解这一问题困难所在。基于此, 笔者对现有的PSO算法、QPSO算法和CQPSO算法进行分析, 并对算法进行修正, 以便更好地适用于复合Bessel多峰函数有限极值点的求解。为了保证上述算法的可靠性, 求取复合Bessel函数零点的同时也采用了逐步搜索二分法进行精细划分求解零点, 将二分法与上述算法求解结果进行对比, 以保证不遗漏零点。

3.1 修正普通粒子群算法

在PSO算法^[8-10]中, 群体内每个个体被称为一个粒子, 表示待求解问题的一个潜在解。粒子i(i=1, 2, …, N)可由2个向量表示, 即速度向量v_{i, t}=[v_{i, t}¹ v_{i, t}² … v_{i, t}^D]和位置向量X_{i, t}=[x_{i, t}¹ x_{i, t}² … x_{i, t}^D], 其中D为待求解问题的维数, t为当前的迭代次数。在每一次迭代计算中, 令向量P_{i, t}=[p_{i, t}¹ p_{i, t}² … p_{i, t}^D]为粒子i到目前搜索到的最优位置, 向量P_{g, t}=[p_{g, t}¹ p_{g, t}² … p_{g, t}^D]为整个群体搜索到的最优位置, 但在PSO算法中, 每一代粒子在进化过程中都以群体最优位置P_{g, t}和个体历史最优位置P_{i, t}加权平均位置为吸引点, 这对于寻求全局最优解很有利。对于函数F(x, r_De)的分析可知, 该多峰函数极小值点全部为0, 若以群体最优位置P_{g, t}和个体历史最优位置P_{i, t}加权位置为吸引点, 则只能求解到部分极值点, 为了增加搜寻有限解空间解的概率, 故将PSO算法中粒子i更新的速度与位置可改变如下

$ x^j_{i,t + 1} = x^j_{i,t} + \omega v^j_{i,t} + {c_1}r^j_1\left( {p^j_{i,t} - x^j_{i,t}} \right), $

(5)

其中: ω为惯性权重; c₁为加速系数; r₁^j为[0, 1]之间满足均匀分布的随机数。通过式(5)可以看出, 每个粒子都以自身历史最优位置为吸引点。

3.2 修正量子粒子群算法

基于对PSO算法中粒子运动轨迹的分析, QPSO算法假定粒子群系统是一个满足量子力学基本假设的粒子系统^[18], 粒子i在第j维上是在以粒子群体最优位置P_{g, t}和个体历史最优位置P_{i, t}加权位置为吸引点, 所有粒子都以吸引点为中心的δ势阱中运动, 具有量子基本行为特征, 则其状态可由波函数ψ描述。为得到粒子的位置, 需要将粒子状态由量子态塌缩到经典态, 采用Monte Carlo随机模拟方法测量粒子的位置, 从而得到粒子的位置更新方程为

$ x^j_{i,t + 1} = \beta ^j_{i,t}p^j_{i,t} + \left( {1 - \beta ^j_{i,t}} \right)p^j_{g,t} \pm a|c^j_t - x^j_{i,t}|\ln\left( { - \frac{1}{{u^j_{i,t}}}} \right), $

(6)

式中:β_{i, t}^j, u_{i, t}^j为(0, 1)之间满足均匀分布的随机数; a为收缩扩张系数; 向量c_t^j=[c_t¹ c_t² … c_t^D]为平均最优位置, 即所有粒子的个体历史最优位置的平均值。

为了增加每个粒子最优解的概率, 减小全局最优解对其他粒子的影响, 更好地寻求有限解空间复合Bessel函数的极值点, 可以将式(6)改变如下:

$ x^j_{i,t + 1} = \beta ^j_{i,t}p^j_{i,t} + \left( {1 - \beta ^j_{i,t}} \right)x^j_{i,t} \pm a|p^j_{i,t} - x^j_{i,t}|\ln\left( { - \frac{1}{{u^j_{i,t}}}} \right)。$

(7)

3.3 修正带交叉算子量子粒子群算法

在修正QPSO算法^[18-19]中, 每个粒子将个体历史最优位置P_{i, t}与当前粒子位置X_{i, t}的加权位置作为自身的吸引点。这种计算方法以粒子运动轨迹的分析结果^[7]为指导, 具有计算简单的优点, 但也存在2个明显的问题: (1)每个粒子除了学习自身的经验外, 只以个体的历史最优位置为引导, 导致计算收敛速度慢, 降低算法求解计算效率; (2)根据式(7)得到的吸引点将限定在以个体历史最优位置P_{i, t}和目前个体位置X_{i, t}为对顶点的D维超矩形内, 随着算法演化过程的进行, 每个粒子吸引点的可能分布空间逐渐减小, 最终导致算法在解空间搜索的能力不足。为了避免吸引点只能限定在特定区域的问题, 可将式(7)改写为

$ x^j_{i,t} = {\beta _{i,t}}p^j_{i,t} + \left( {1 - {\beta _{i,t}}} \right)x^j_{i,t} + \frac{{p^j_{b,t} - p^j_{c,t}}}{2} + a\left| {\frac{{p^j_{b,t} - p^j_{c,t}}}{2}} \right|\ln\left( { - \frac{1}{{u^j_{i,t}}}} \right), $

(8)

式中:下标b, c为群体中随机选取的2个粒子编号, 且b≠c≠i。x_{i, t}^j需要进行交叉操作, 具体的交叉操作可参见文献[18], 通过这种关于粒子吸引点的计算方法, 可以提高算法求解多峰函数优化问题的性能。这是因为当大部分粒子集聚在一起时, 分布在其他搜索空间中的粒子有机会帮助部分停滞粒子逃离进入到其他可行解区域。

4 复合Bessel函数零点分布规律 4.1 复合Bessel函数零点计算结果

在PSO算法中, 惯性权ω∈[0.4, 0.95], 采用线性递减策略, 惯性权随着迭代次数的增加而减小。c₁加速系数取1。文献[19]指出, 收缩扩张系数a的取值方法与取值范围会对粒子的行为产生影响。a取值越大, 粒子找到新知识的可能性越大、创造力越强; a取值越小, 粒子的局部挖掘开发能力越强、精细搜索效果越好。修正QPSO和修正CQPSO算法中的收缩扩张系数a, 通过多次试验, 结合文献[15], 采用“Z”型递减策略, 在早期, 增加粒子探索范围, 在后期递减, 增加粒子的收敛速度, 具体收缩扩张系数a计算公式如下

$ a(t)=\frac{1.58}{1+\exp(0.02t)}, $

(9)

其中:t为迭代次数。

选取最大迭代次数为2 000, 种群粒子个数统一设置为1 000。为了减小统计误差, 3种算法在计算机上独立运行50次, 取50次运行的平均结果进行分析与比较。试验的测试环境设置如下:操作系统为Windows 7, Intel(R) Xeon(R) CPU E3-1230 v3, 3.30 GHz 8GB计算机硬件环境条件下, 利用Matlab 2013b软件编写3种算法程序, 分别运行独立运行50次。

鉴于篇幅有限, 这里只给出了r_De=10时, 搜索范围为[0.000 000 1, 3.5], 通过二分法, 修正PSO法, 修正QPSO法和修正CQPSO法计算得到前10个零点。通过表 1可以看出, 将修正PSO、QPSO和CQPSO计算结果与二分法进行对比, 多次独立运行修正这3种算法，都能找到搜索范围内所有实根，在相同环境中，修正PSO和修正CQPSO算法计算的精度较修正QPSO算法高, 且修正CQPSO算法独立运行50次, 搜索到零点的概率为100%, 而修正PSO和修正QPSO算法对方程前3个零点有很小概率搜索不到, 这说明3种算法通过修正后进行多次独立运行都适用于复合Bessel多峰函数零点的求解。

图 2为3种算法求解复合Bessel函数零点求解的收敛曲线, 图 2为独立计算50次的绝对误差与迭代次数关系, 从图 2中可以看出, 修正CQPSO和修正PSO算法收敛稳定, 而修正QPSO收敛相对不太稳定。为便于3种算法收敛速度的对比, 选取迭代2 000次, 可以明显看出, 修正CQPSO算法比修正QPSO和修正PSO算法收敛速度快得多, 在迭代次数相同的情况下, 前者比后者求解零点精度高, 而修正QPSO算法收敛速度慢, 且精度低, 为此取修正CQPSO和修正PSO算法中精度较高者作为函数零点。

4.2 复合Bessel函数零点分布规律

对不同参数的复合Bessel函数零点采用二分法、修正PSO和修正CQPSO算法进行求解, 将不同参数求解得到的零点与零点次序绘制在双对数坐标系中, 见图 3。除前3个零点之外, 对于不同参数的复合Bessel函数零点(第4个零点以后)与零点次序呈直线关系。从图 3可以看出, 当r_De越大时, 函数零点越小, 即在相同求解区间内零点分布个数越多, 这与图 1中零点分布规律一致。

为了确定复合Bessel函数零点与零点次序线性关系, 结合复合Bessel方程形式, 不妨假设直线方程为

$ \lg\left( {{\alpha _i}} \right) = k\lg\left( i \right) + m, $

(10)

其中:k为直线的斜率, m为直线截距, i(i=4, 5, …)为零点次序, α_i(i=4, 5, …)为方程的零点。

对求解得到不同参数情况下函数零点与零点次序进行直线方程进行拟合, 拟合结果见表 2, 除去前3个零点之外, 对于同一参数的复合Bessel函数零点与零点次序进行直线方程拟合, 拟合的相关系数基本都接近于1, 相关程度高。将拟合系数带入直线方程计算得到函数拟合零点, 把拟合零点带入到复合Bessel方程计算的函数值与真值之差都小于10^-12。通过拟合零点与求解零点计算相对误差, 由表 2可以看出, 本研究计算的结果与Hurst计算的结果一致, 计算的相对误差都小于0.5%, 能够满足工程计算的需求。在具体工程应用中, 只要得到面前几个零点, 可通过拟合直线方程来计算后面的零点, 简化复杂Bessel函数的计算。

5 结论及建议

(1) 复合Bessel函数f(x, r_De)是以r_De为参数的多峰函数, 函数零点有无穷多个, 当参数r_De取值较小时, 零点值较大, 且分布相对较稀疏; 而当参数r_De取值较大时, 零点值较小, 且分布密集。在零点求取过程中, 若参数r_De较小, 可适当放大搜索范围; 反之, 则要缩小搜索范围。

(2) 对PSO算法和QPSO算法进行个体历史最优位置和群体最优位置加权为吸引点的计算方法进行修正, 同时对QPSO算法进行分析基础之上对CQPSO算法进行修正, 将这3种修正算法求解的结果与二分法逐步搜索算法进行对比, 说明3种算法在求解复合Bessel函数零点时是可行的。

(3) 通过修正CQPSO算法与其他两种算对比, 修正CQPSO算法收敛速度快, 且计算精度高, 说明修正CQPSO算法中的交叉算子操作在优化多峰函数时能够取得良好的效果。

(4) 复合Bessel函数除去前3个零点, 后续零点与零点次序在双对数坐标系中满足直线关系, 对同一参数r_De拟合结果表明:零点与零点次序的线性相关程度高, 拟合零点与求解零点相对误差不超过0.5%, 因此可通过给出3种优化算法获取任意参数r_De前面有限个零点, 然后通过直线方程拟合可以获取函数其他任意零点, 在工程计算避免复杂的Bessel函数的计算。