CN103886221A - 一种变质量破碎岩石渗透性参量计算方法 - Google Patents

Abstract

本发明涉及一种变质量破碎岩石渗透性参量计算方法,该方法设计一种遗传算法，对渗透率、非Darcy流β因子、加速度系数三者之间幂指数函数关系和中渗透率的参考值k_r、非Darcy流β因子的参考值β_r、加速度系数的参考值c_ar、幂指数n_β和n_c进行优化，可以克服试验者经验不足的困境，也不需要参考先前的试验结果，即可获取合理的决策变量；可以输出初始种群及各代新种群的位串、k_r、β_r、c_arn_β和n_c、渗流速度时间曲线、最佳个体的渗透率、非Darcy流β、加速度系数和适应度；根据压力梯度和渗流速度时间序列计算各采样时刻的渗透率、非Darcy流β因子和加速度系数，为分析质量流失引起的破碎岩石渗透性参量的变化规律提供条件。

Description

一种变质量破碎岩石渗透性参量计算方法

技术领域

本发明涉及一种岩石力学性质测试的计算方法，具体涉及一种变质量破碎岩石渗透性参量计算方法。

背景技术

在破碎岩石渗透试验过程中，由于水的溶蚀、磨蚀和冲蚀作用，细小颗粒从破碎岩石的表面分离出来，并在破碎岩石的孔隙中迁移。细小颗粒的迁移造成破碎岩石的质量流失，故破碎岩石的孔隙度和渗透性参量（渗透率、非Darcy流β因子、加速度系数）随时间变化。借助于实验室试验，分析破碎岩石在渗透过程中质量迁移与流失规律、孔隙度和渗透性参量时变规律，对于破碎岩体渗流系统的数字仿真、渗流失稳的机理分析具有基础性和支撑性作用。根据压力梯度和渗流速度时间序列计算各采样时刻的渗透率、非Darcy流β因子和加速度系数，是破碎岩体渗流系统非线性动力学行为的重要环节。

目前，变质量破碎岩石渗透性参量参照恒定质量破碎岩石的方法进行，对于质量慢变的破碎岩石来说，精度一般可以达到工程设计的要求，但当质量流失较大时，计算的渗透性参量误差较大。

目前，计算变质量破碎岩石在采样时刻t_i＝iτ(i＝1,2,…,N;τ为采样周期)的渗透率、非Darcy流β因子和加速度系数是建立在渗透率、非Darcy流β因子、加速度系数三者之间满足幂指数函数关系的基础上的，即将恒定质量破碎岩石在确定的孔隙度和确定的孔隙度下的渗透率、非Darcy流β因子、加速度系数作为参考值，而是时变孔隙度和时变压力梯度下的渗透率、非Darcy流β因子、加速度系数三者之间满足如下的幂指数关系

$\mathrm{\β}={\mathrm{\β}}_r{\left(\frac{k}{k_r}\right)}^{n_{\mathrm{\β}}}---(a-1)$

$c_a=c_{\mathrm{ar}}{\left(\frac{k}{k_r}\right)}^{n_c}---(a-2)$

其中，幂指数n_β和n_c也参照恒定质量破碎岩石的渗透试验结果确定。具体的做法分为五步：第一步，构造时间序列，包括渗流速度时间序列V_i(i＝0,1,2,L,N)、压力梯度时间序列

和当地加速度

的时间序列a_i(i＝0,1,2,L,N)；

第二步，建立采样时刻渗透率的代数方程，即将动量守恒方程与式（a-1）、式（a-2）结合，得到：

$\mathrm{\ρ}{c}_a^r{\left(\frac{k_i}{k_r}\right)}^{m_c}{a}_i=-G_i-\frac{\mathrm{\μ}}{k_i}{V}_i-\mathrm{\ρ}{\mathrm{\β}}_r{\left(\frac{k_i}{k_r}\right)}^{m_{\mathrm{\β}}}{V}_i^2,i=\mathrm{0,1,2},...N---(a-3)$

利用Newton切线法代数方程（a-3）的根，得到采样时刻的渗透率，再利用式（a-2）和式（a-3）分别计算出采样时刻的非Darcy流β因子和加速度系数；

第三步，构造动量守恒方程的外部函数，利用变步长的四阶Runge-Kutta法求方程动量守恒方程：

$\mathrm{\ρ}{c}_a\frac{\∂V}{\∂t}=-\frac{\∂p}{\∂x}-\frac{\mathrm{\μ}}{k}V-\mathrm{\ρ\β}{V}^2---(a-4)$

的数值解，需要构造压力梯度、渗透率、非Darcy流β因子和加速度系数非采样时刻（如t＝2.2τ）的值，在变步长的四阶Runge-Kutta法子程序中，外部函数需要调用全区间三节点插值子程序，分别对压力梯度、渗透率、非Darcy流β因子和加速度系数进行全区间三节点插值。例如，时刻t∈(t_k,t_k+2),(k＝0,1,2,L,N)的渗透率，可按如下的插值公式计算：

$\begin{array}{c}k\left(t\right)=k_k\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}\\ +k_{k+1}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}\\ +k_{k+2}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}\end{array};---(a-5)$

第四步，以τ为步长，利用变步长的四阶Runge-Kutta法求方程（a-1）的数值解，得到渗流速度计算值时间序列。

第五步，计算渗流速度的数值解

与试验数据V_i(i＝0,1,2,L,N)之间的差E_rr，

$E_{\mathrm{rr}}=\frac{1}{N}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(1-\frac{V_i^{*}}{V_i})}^2---(a-6)$

第六步，如果E_rr≤ε（ε为某一预设值），则可认为参数合理，并将在这组参数下计算得到的渗透率、非Darcy流β因子和加速度系数序列作为渗透性参量的测试结果。如果E_rr>ε，则需要对参数进行调整，再重复第一至第四步。

目前的计算方法不足之处在于渗透率、非Darcy流β因子、加速度系数的参考值、幂指数n_β和n_c的选择完全依赖于试验者的经验和先前的试验数据。由于不同矿区的破碎岩石的渗透性差异很大，参考先前的试验数据确定决策变量k_r、β_r、c_arn_β和n_c的取值，可能出现较大误差。试验者的经验积累也是漫长的过程，只有短时间试验经历的试验者难以找到合理的决策变量。

发明内容

为了解决上述问题，本发明提出了一种变质量破碎岩石渗透性参量计算方法，该方法设计一种遗传算法，对渗透率、非Darcy流β因子、加速度系数三者之间幂指数函数关系

和

中渗透率的参考值k_r、非Darcy流β因子的参考值β_r、加速度系数的参考值c_ar、幂指数n_β和n_c进行优化，可以克服试验者经验不足的困境，也不需要参考先前的试验结果，即可获取合理的决策变量；可以输出初始种群及各代新种群的位串、k_r、β_r、c_arn_β和n_c、渗流速度时间曲线、最佳个体的渗透率、非Darcy流β、加速度系数和适应度；根据压力梯度和渗流速度时间序列计算各采样时刻的渗透率、非Darcy流β因子和加速度系数，为分析质量流失引起的破碎岩石渗透性参量的变化规律提供条件。

本发明是通过以下技术方案实现的：一种变质量破碎岩石渗透性参量计算方法，是变质量破碎岩石在渗透试验过程中各采样时刻的渗透性参量，包括渗透率k、非Darcy流β因子β、加速度系数c_a的计算方法，具体操作步骤为：

步骤1，假设渗透率k、非Darcy流β因子β、加速度系数c_a之间存在幂指数关系，即，

$\mathrm{\β}={\mathrm{\β}}_r{\left(\frac{k}{k_r}\right)}^{n_{\mathrm{\β}}},---(1-1)$

$c_a=c_{\mathrm{ar}}{\left(\frac{k}{k_r}\right)}^{n_c};---(1-2)$

步骤2，构造压力梯度、渗流速度和渗流速度当地导数的时间序列，采样周期为τ，长度为N的压力梯度G_p、渗流速度V和渗流速度当地导数的时间序列分别记为：V_i(i＝0,1,2,L,N)，a_i(i＝0,1,2,L,N)；

步骤3，对渗透率、非Darcy流β因子、加速度系数的参考值k_r、β_r、c_ar、幂指数n_β和n_c进行编码及解码；

步骤3.1，参照恒定质量破碎岩石渗透试验的结果，确定渗透率、非Darcy流β因子、加速度系数的参考值k_r、β_r、c_ar、幂指数n_β和n_c的可能取值范围；

步骤3.2，确定基因长度，将决策变量渗透率、非Darcy流β因子、加速度系数的参考值k_r、β_r、c_ar、幂指数n_β和n_c，转换为长度分别为l₁、l₂、l₃、l₄和l₅，由字符0和1组成的位串

和

；

步骤3.3，由长度为m＝l₁+l₂+l₃+l₄+l₅的二进制位串I₁I₂…I_m构成遗传算法的个体基因型；

步骤3.4，产生初始种群，并对种群中每一个体进行解码操作，计算出个体的k_r、β_r、c_ar、n_β和n_c的值；

步骤3.4.1，确定初始种群规模k_group；

步骤3.4.2，产生随机种子数；

步骤3.4.3，产生k_group个长度为m的二进制位串，得到初始种群

$\mathrm{initial}\_\mathrm{pop}=\{I_1^i{I}_2^i...I_m^i|i=\mathrm{1,2},...,k_{\mathrm{group}}\};$

步骤3.4.4，对种群中每一个体进行解码操作，得到个体基因的表现型，计算出个体的k_r、β_r、c_ar、n_β和n_c的值；

步骤4，获取种群中每一个体对应的渗透率、非Darcy流β因子和加速度系数时间序列，构造动量守恒方程的外部函数，将渗流速度计算值与试验数据进行对比，计算种群中个体的适应度；

步骤4.1，利用Newton切线法求采样时刻的渗透率、非Darcy流β因子和加速度系数，得到每一个体对应的渗透率、非Darcy流β因子和加速度系数时间序列，

${\mathrm{\β}}^{t_i}(i=\mathrm{1,2},...,N),c_a^{\mathrm{ti}}=(i=\mathrm{1,2},...,N);$

步骤4.2，分别对压力梯度、渗透率、非Darcy流β因子和加速度系数进行全区间三节点插值构造动量守恒方程的外部函数；

步骤4.3，以τ为步长，利用变步长的四阶Runge-Kutta法求动量守恒方程的数值解，得到采样时刻的渗流速度计算值时间序列

步骤4.4，计算渗流速度的计算值

与试验数据V_i(i＝0,1,2,L,N)之间的误差E_rr；

步骤4.5，构造适应度函数并以误差E_rr计算种群中个体的适应度；

步骤5，通过对种群进行选择、交叉、变异操作获得新的种群，对新种群中的个体的适应度进行排序以得到最佳个体，并最佳个体进行解码；利用Newton切线法计算，获取最佳个体的渗透率、非Darcy流β因子和加速度系数，具体操作步骤为：

步骤5.1，选择操作，从初始种群中选择出具有交配权的k_cop个体；

步骤5.2，交叉操作，对每一基因位串随机产生交叉位，并按某一交叉概率p_c对每一对交叉个体（夫妇）进行交叉操作，得到个体集合；

步骤5.3，变异操作，对个体集合中每一基因位串随机产生变异位，并按某一交叉概率p_c对每一个体进行变异操作，得到新一代种群；

步骤5.4，对新种群中的个体进行解码操作，计算每一个体的适应度；

步骤5.5，对新种群中的每一个体适应度进行排序；

步骤5.6，设置停止繁殖的条件；

步骤5.7，对最佳个体（适应度最大的个体）进行解码，将最佳个体的基因型转化为表现型，得到k_r、β_r、c_arn_β和n_c的估计值；

步骤5.8，利用Newton切线法计算出最佳个体各采样时刻的渗透率、非Darcy流β因子和加速度系数。

进一步地，步骤3.1中，渗透率、非Darcy流β因子、加速度系数的参考值k_r、β_r、c_ar、幂指数n_β和n_c的可能取值范围是：

k_r∈[k_r1,k_r2]，    （3-1）

β_r∈[β_r1,β_r2]，    （3-2）

c_ar∈[c_ar1,c_ar2]，    （3-3）

n_β∈[n_β1,n_β2]，    （3-4）

n_c∈[n_c1,n_c2]。    （3-5）

进一步地，步骤3.3中所述遗传算法的个体基因型，相应的表现型为：

$k_r^j=k_{r1}{\left[\exp \frac{\ln \frac{k_{r2}}{k_{r1}}}{2^{l_1}-1}\right]}^j,j=\underset{i=1}{\overset{l_1}{\mathrm{\Σ}}}{2}^i{I}_{1i},j=\mathrm{0,1},...,2^{l_1}-1---(3-6)$

${\mathrm{\β}}_r^j={\mathrm{\β}}_{r1}{\left[\exp \frac{\ln \frac{{\mathrm{\β}}_{r2}}{{\mathrm{\β}}_{r1}}}{2^{l_2}-1}\right]}^j,j=\underset{i=1}{\overset{l_2}{\mathrm{\Σ}}}{2}^i{I}_{2i},j=\mathrm{0,1},...,2^{l_2}-1---(3-7)$

$c_{\mathrm{ar}}^j=c_{\mathrm{ar}1}{\left[\exp \frac{\ln \frac{c_{\mathrm{ar}2}}{c_{\mathrm{ar}1}}}{2^{l_3}-1}\right]}^j,j=\underset{i=1}{\overset{l_3}{\mathrm{\Σ}}}{2}^i{I}_{3i},j=\mathrm{0,1},...,2^{l_3}-1---(3-8)$

$n_b^j=n_{b1}[1+\frac{j}{2^{l_4}-1}\frac{n_{b2}-n_{b1}}{n_{b1}}],j=\underset{i=1}{\overset{l_4}{\mathrm{\Σ}}}{2}^i{I}_{4i},j=\mathrm{0,1},...,2^{l_4}-1---(3-9)$

$n_c^j=n_{c1}[1+\frac{j}{2^{l_5}-1}\frac{n_{c2}-n_{c1}}{n_{c1}}],j=\underset{i=1}{\overset{l_5}{\mathrm{\Σ}}}{2}^i{I}_{5i},j=\mathrm{0,1},...,2^{l_5}-1---(3-10)$

进一步地，步骤3.4.4中所述的个体基因的表现型为：

$k_r^j=k_{r1}\exp \left(\frac{j}{2^{l_1}-1}\ln \frac{k_{r2}}{k_{r1}}\right),j=\underset{i=1}{\overset{l_1}{\mathrm{\Σ}}}{2}^i{I}_{1i},j=\mathrm{0,1},...,2^{l_1}-1---(3-11)$

${\mathrm{\β}}_r^j={\mathrm{\β}}_{r1}\exp \left(\frac{j}{2^{l_2}-1}\ln \frac{{\mathrm{\β}}_{r2}}{{\mathrm{\β}}_{r1}}\right),j=\underset{i=1}{\overset{l_2}{\mathrm{\Σ}}}{2}^i{I}_{2i},j=\mathrm{0,1},...,2^{l_2}-1---(3-12)$

$c_{\mathrm{ar}}^j=c_{\mathrm{ar}1}\exp \frac{j}{2^{l_3}-1}\ln \frac{c_{\mathrm{ar}2}}{c_{\mathrm{ar}1}},j=\underset{i=1}{\overset{l_3}{\mathrm{\Σ}}}{2}^i{I}_{3i},j=\mathrm{0,1},...,2^{l_3}-1---(3-13)$

$n_b^j=n_{b1}[1+\frac{j}{2^{l_4}-1}\frac{n_{b2}-n_{b1}}{n_{b1}}],j=\underset{i=1}{\overset{l_4}{\mathrm{\Σ}}}{2}^i{I}_{4i},j=\mathrm{0,1},...,2^{l_4}-1---(3-14)$

$n_c^j=n_{c1}[1+\frac{j}{2^{l_5}-1}\frac{n_{c2}-n_{c1}}{n_{c1}}],j=\underset{i=1}{\overset{l_5}{\mathrm{\Σ}}}{2}^i{I}_{5i},j=\mathrm{0,1},...,2^{l_5}-1---(3-15)$

进一步地，步骤4.2中所述的动量守恒方程的外部函数为：

$\frac{\∂V}{\∂t}=-\frac{1}{{\mathrm{\ρc}}_a}(G_p+\frac{\mathrm{\μ}}{k}V+\mathrm{\ρ\β}{V}^2)---(4-1)$

其中，ρ为水的质量密度，β为水的动力粘度；其中，t∈[t_k,t_k+2]时刻的渗透率、非Darcy流β因子、加速度系数和压力梯度按全区间不等距三节点Lagrange插值计算如下：

$k=k^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+k^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+k^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---(4-2)$

$\mathrm{\β}={\mathrm{\β}}^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+{\mathrm{\β}}^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+{\mathrm{\β}}^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---(4-3)$

$c_a=c_a^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+c_a^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+c_k^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---(4-4)$

$G_p=G_p^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+G_p^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+G_p^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}.---(4-5)$

进一步地，步骤4.4中所述的E_rr的计算公式为： $E_{\mathrm{rr}}=\frac{1}{N}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(1-\frac{V_i^{*}}{V_i})}^2.$ $(4-6)$

进一步地，步骤4.5中所述的计算种群中个体的适应度的公式为：

$\mathrm{fitn}\left(i\right)=\frac{1}{{\mathrm{err}}^{1/4}},i=\mathrm{1,2},...,k_{\mathrm{group}}.---(4-7)$

进一步地，步骤5.1中所述的具有交配权的k_cop个体的基因型为：

$P_{\mathrm{cop}}=\left\{\mathrm{chromosome}\left(i_1\right)\right|\mathrm{chromosome}\left(i_1\right)=I_1^{i_1}{I}_2^{i_1}...I_m^{i_1},i=\mathrm{1,2},...,k_{\mathrm{cop}}\}.---(5-1)$

进一步地，经步骤5.2进行交叉运算后得到个体集合

；经步骤5.3对

中每一个体进行变异运算得到新一代种群： $\mathrm{NEW}\_\mathrm{pop}=\{I_1^i{I}_2^i...I_m^i|i=\mathrm{1,2},...,k_{\mathrm{group}}\}.---(5-2)$

进一步地，步骤5.4中所述的计算新一代种群中每一个体的适应度fitn(i)，i＝1,2,…,k_group。

进一步地，步骤5.5中，如果新一代群体中个体适应度的最大值fitn_max大于等于预先设定的数值s，即fitn_max≥s则停止繁殖；如果fitn_max＜s，则继续进行选择、交叉、变异运算，直到fitn_max≥s得以满足或繁殖代数等于某一事先设定的值T，则停止繁殖。

以上所述步骤5.4所用的公式同步骤4.5所用的公式相同；步骤5.7所用的公式同步骤3.4.4所用的公式相同。

本发明的有益效果是：一种变质量破碎岩石渗透性参量计算方法，具有以下优势，（1）利用遗传算法对参数k_r、β_r、c_arn_β和n_c进行优化，可以克服试验者经验不足的困境，也不需要参考先前的试验结果，即可获取合理的决策变量；（2）试验者只需将渗透过程中压力（差）和流量进行周期采样，产生一个包含三列数据（第一列为时间，第二列为压力/压力差，第三列为流量）的文本文件或Data文件，便可在程序中修改某些参数进行参数优化和渗透率计算，这些参数包括液体的质量密度、动力粘度、岩样直径、厚度、质量、采样周期、数据长度、交叉概率、变异概率、初始种群规模、交叉个体对数、交叉位、变异位、基因长度等；（3）可以输出初始种群及各代新种群的位串、k_r、β_r、c_arn_β和n_c、渗流速度-时间曲线、最佳个体的渗透率、非Darcy流β、加速度系数和适应度；（4）根据压力梯度和渗流速度时间序列计算各采样时刻的渗透率、非Darcy流β因子和加速度系数，其精确度高，计算出的渗透性参量误差较小，为分析质量流失引起的破碎岩石渗透性参量的变化规律提供了有利条件。

附图说明

下面结合附图及实施例对本发明作进一步说明。

附图1为压力梯度和渗流速度的时间序列图。

附图2是渗透率的时间序列图。

附图3是非Darcy流β因子的时间序列图。

附图4是加速度系数的时间序列图。

附图5渗流速度的实测值与计算值的时间序列图。

具体实施方式

利用本发明进行变质量破碎岩石渗透性参量计算，其具体操作如下：

1.试验数据处理：

将试验数据改为三列，第一列为时间，单位为s；第二列为压力（差），单位为MPa；第三列为流量，单位是l/h。

2.输入有关参数：

输入有关参数，包括以下几类：

（1）液体的质量密度ρ、动力粘度μ；

（2）岩样半径r_s、自然堆放高度h₀、压缩后的高度h、质量；

（3）采样周期τ、数据长度N；

（4）交叉概率、变异概率、初始种群规模、交叉个体对数、交叉位、变异位、基因长度等；

（5）描述停止繁殖的条件的参数，如最大繁殖代数T、最小适应度s；

（6）Newton切线法的最大允许误差ε₁、Runge-Kutta法的最大允许误差ε；

（7）描述k_r、β_r、c_arn_β和n_c取值范围的参数；

本例中，这些参数分别为：

（1）液体的质量密度ρ＝1000(kg/m³)、动力粘度μ＝1.01×10^-3(Pa·s)；

（2）岩样半径r_s＝0.25×10^-2(m)、自然堆放高度h₀＝1.25×10^-1(m)、压缩后的高度h＝1.00×10^-1(m)、质量m＝1800(kg)；

（3）采样周期τ＝1.0(s)、数据长度N＝301；

（4）交叉概率p_c＝0.8、变异概率p_m＝0.2、初始种群规模k_group＝40、交叉个体对数k_cop＝10，基因k_r、β_r、c_arn_β和n_c的长度分别为5、5、5、4和4，交叉位和变异位随机产生。

（5）最大繁殖代数T＝31、最小适应度s＝2000；

（6）Newton切线法的最大允许误差ε₁＝1.0×10^-9、Runge-Kutta法的最大允许误差ε＝1.0×10^-4；

（7）描述k_r、β_r、c_arn_β和n_c取值范围的参数：

k_r1＝0.5×10^-12(m²),k_r2＝0.7×10^-12(m²)，

β_r1＝0.24×10⁺¹⁰(m^-1),k_r2＝0.26×10⁺¹⁰(m^-1)，

c_ar1＝0.18×10⁺¹⁰,c_ar2＝0.22×10⁺¹⁰，

n_β1＝-0.21×10⁺⁰¹,n_β2＝-0.17×10⁺⁰¹，

n_c1＝-0.15×10⁺⁰¹,n_c2＝-0.11×10⁺⁰¹。

3.压力梯度和渗流速度时间序列的产生：

读取试验中压力（差）和流量采样数据文件，经过简单运算，得到压力梯度和渗流速度的时间序列，见图1。

4.初始种群的产生

利用Fortran程序随机生成规模为40的初始种群，见表1。

表1初始种群

5.繁殖130909000

利用随机遍历法从初始种群中选择出具有交配权的20个体，对20中个体进行随机配对；对每一基因随机产生交叉位，并按p_c＝0.8的交叉概率进行交叉操作。对交叉操作后的群体中每一个体五个基因随机产生变异位，并按p_m＝0.2的变异概率进行变异操作。对变异操作后的群体中个体进行解码操作，计算相应的适应度；对适应度进行排序，得到适应度最大的个体序号k_best。如果fitn(k_best)＜s且繁殖代数i_g＜T，则继续繁殖，直到fitn(k_best)＞s或繁殖代数i_g＝T，则停止繁殖；对编号为k_best的个体进行解码，得到最优基因

和

经过两代繁殖，便可得到fitn(k_best)＞s＝2000的个体（基因）。第一代和第二代个体基因及适应度分别见表2和表3。可见，第二代群体中第8个体的适应度为2120，大于2000。第8个体的基因即为最优基因，且

${\mathrm{\β}}_r^{\mathrm{best}}=2.42\×{10}^{+09}\left(m^{-1}\right),$ $c_{\mathrm{ar}}^{\mathrm{best}}=1.85\×{10}^{+09},$ $n_{\mathrm{\β}}^{\mathrm{best}}=-2.02$ 和 $n_c^{\mathrm{best}}=-1.50.$

表2第一代个体基因及适应度

序号

kr

βr

car

nβ

nc

适应度

1

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

2

5.22E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.90E+03

3

5.00E-13

2.42E+09

1.85E+09

-1.97E+00

-1.50E+00

1.84E+03

4

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

5

5.00E-13

2.43E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

6

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

7

5.00E-13

2.42E+09

2.05E+09

-2.02E+00

-1.50E+00

1.94E+03

8

5.00E-13

2.53E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

9

5.00E-13

2.48E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

10

5.00E-13

2.42E+09

1.85E+09

-1.97E+00

-1.50E+00

1.84E+03

11

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

12

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

13

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

14

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

15

5.22E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.90E+03

16

5.00E-13

2.42E+09

1.85E+09

-1.81E+00

-1.50E+00

1.84E+03

17

5.00E-13

2.42E+09

2.05E+09

-2.02E+00

-1.50E+00

1.94E+03

18

5.34E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.93E+03

19

5.00E-13

2.42E+09

1.85E+09

-1.99E+00

-1.50E+00

1.84E+03

20

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

表3第二代个体基因及适应度

序号

kr

βr

car

nβ

nc

适应度

1

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

2

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

3

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

4

5.00E-13

2.42E+09

1.85E+09

-1.81E+00

-1.50E+00

1.84E+03

5

5.05E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.86E+03

6

5.11E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.87E+03

7

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

8

5.95E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

2.12E+03

9

5.28E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.91E+03

10

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

11

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

12

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

13

5.00E-13

2.53E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

14

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

15

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

16

5.22E-13

2.42E+09

2.05E+09

-2.02E+00

-1.50E+00

2.01E+03

17

5.00E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

18

5.00E-13

2.43E+09

1.85E+09

-2.02E+00

-1.50E+00

1.84E+03

19

5.00E-13

2.42E+09

1.85E+09

-1.97E+00

-1.50E+00

1.84E+03

20

5.34E-13

2.42E+09

1.85E+09

-2.02E+00

-1.50E+00

1.93E+03

6.渗透性参量计算

针对最优的基因

和

利用Newton切线法代数方程求解采样时刻渗透率满足的代数方程：

$\mathrm{\ρ}{c}_{\mathrm{ar}}^{\mathrm{best}}{\left(\frac{k_i}{k_r}\right)}^{n_c^{\mathrm{best}}}{a}_i=-G_i-\frac{\mathrm{\μ}}{k_i}{V}_i-\mathrm{\ρ}{\mathrm{\β}}_r^{\mathrm{best}}{\left(\frac{k_i}{k_r^{\mathrm{best}}}\right)}^{n_{\mathrm{\β}}^{\mathrm{best}}}{V}_i^2,i=\mathrm{0,1,2},...N$

得到采样时刻的渗透率，再分别计算出采样时刻的非Darcy流β因子和加速度系数。这样得到渗透率、非Darcy流β因子和加速度系数的时间序列，见图2～4。

7.渗流速度和计算值与测试值的比较

基于最优基因 $k_r^{\mathrm{best}}=5.95\×{10}^{-13}\left(m^2\right),$ ${\mathrm{\β}}_r^{\mathrm{best}}=2.42\×{10}^{+09}\left(m^{-1}\right),$ $c_{\mathrm{ar}}^{\mathrm{best}}=1.85\×{10}^{+09},$ 和

求动量守恒方程的数值解，得到渗流速度计算值时间序列

根据fitn(k_best)反算出渗流速度计算值与实测值的误差，有

$\mathrm{err}={\left[\frac{1}{\mathrm{fitn}\left(k_{\mathrm{best}}\right)}\right]}^4={\left(\frac{1}{2120}\right)}^4=4.95\×{10}^{-14}$

渗流速度的实测值与计算值的时间序列，见图5。

Claims (10)

1.一种变质量破碎岩石渗透性参量计算方法，是变质量破碎岩石在渗透试验过程中各采样时刻的渗透性参量，包括渗透率k、非Darcy流β因子β、加速度系数c_a的计算方法，其特征在于，具体操作步骤为：

步骤1，假设渗透率k、非Darcy流β因子β、加速度系数c_a之间存在幂指数关系，即，

$\mathrm{\β}={\mathrm{\β}}_r{\left(\frac{k}{k_r}\right)}^{n_{\mathrm{\β}}},---(1-1)$

$c_a=c_{\mathrm{ar}}{\left(\frac{k}{k_r}\right)}^{n_c};---(1-2)$

步骤2，构造压力梯度、渗流速度和渗流速度当地导数的时间序列，采样周期为τ，长度为N的压力梯度G_p、渗流速度V和渗流速度当地导数

的时间序列分别记为：

V_i(i＝0,1,2,L,N)，a_i(i＝0,1,2,L,N)；

步骤3，对渗透率、非Darcy流β因子、加速度系数的参考值k_r、β_r、c_ar、幂指数n_β和n_c进行编码及解码；

步骤3.1，参照恒定质量破碎岩石渗透试验的结果，确定渗透率、非Darcy流β因子、加速度系数的参考值k_r、β_r、c_ar、幂指数n_β和n_c的可能取值范围；

步骤3.2，确定基因长度，将决策变量渗透率、非Darcy流β因子、加速度系数的参考值k_r、β_r、c_ar、幂指数n_β和n_c，转换为长度分别为l₁、l₂、l₃、l₄和l₅，由字符0和1组成的位串

和

；

步骤3.3，由长度为m＝l₁+l₂+l₃+l₄+l₅的二进制位串I₁I₂…I_m构成遗传算法的个体基因型；

步骤3.4，产生初始种群，并对种群中每一个体进行解码操作，计算出个体的k_r、β_r、c_ar、n_β和n_c的值；

步骤3.4.1，确定初始种群规模k_group；

步骤3.4.2，产生随机种子数；

步骤3.4.3，产生k_group个长度为m的二进制位串，得到初始种群

$\mathrm{initial}\_\mathrm{pop}=\{I_1^i{I}_2^i...I_m^i|i=\mathrm{1,2},...,k_{\mathrm{group}}\};$

步骤3.4.4，对种群中每一个体进行解码操作，得到个体基因的表现型，计算出个体的k_r、β_r、c_ar、n_β和n_c的值；

步骤4，获取种群中每一个体对应的渗透率、非Darcy流β因子和加速度系数时间序列，构造动量守恒方程的外部函数，将渗流速度计算值与试验数据进行对比，计算种群中个体的适应度；

步骤4.1，利用Newton切线法求采样时刻的渗透率、非Darcy流β因子和加速度系数，得到每一个体对应的渗透率、非Darcy流β因子和加速度系数时间序列，

${\mathrm{\β}}^{t_i}(i=\mathrm{1,2},...,N),c_a^{\mathrm{ti}}=(i=\mathrm{1,2},...,N);$

步骤4.2，分别对压力梯度、渗透率、非Darcy流β因子和加速度系数进行全区间三节点插值构造动量守恒方程的外部函数；

步骤4.3，以τ为步长，利用变步长的四阶Runge-Kutta法求动量守恒方程的数值解，得到采样时刻的渗流速度计算值时间序列

步骤4.4，计算渗流速度的计算值

与试验数据V_i(i＝0,1,2,L,N)之间的误差E_rr；

步骤4.5，构造适应度函数并以误差E_rr计算种群中个体的适应度；

步骤5，通过对种群进行选择、交叉、变异操作获得新的种群，对新种群中的个体的适应度进行排序以得到最佳个体，并最佳个体进行解码；利用Newton切线法计算，获取最佳个体的渗透率、非Darcy流β因子和加速度系数，具体操作步骤为：

步骤5.1，选择操作，从初始种群中选择出具有交配权的k_cop个体；

步骤5.2，交叉操作，对每一基因位串随机产生交叉位，并按某一交叉概率p_c对每一对交叉个体（夫妇）进行交叉操作，得到个体集合；

步骤5.3，变异操作，对个体集合中每一基因位串随机产生变异位，并按某一交叉概率p_c对每一个体进行变异操作，得到新一代种群；

步骤5.4，对新种群中的个体进行解码操作，计算每一个体的适应度；

步骤5.5，对新种群中的每一个体适应度进行排序；

步骤5.6，设置停止繁殖的条件；

步骤5.7，对最佳个体（适应度最大的个体）进行解码，将最佳个体的基因型转化为表现型，得到k_r、β_r、c_arn_β和n_c的估计值；

步骤5.8，利用Newton切线法计算出最佳个体各采样时刻的渗透率、非Darcy流β因子和加速度系数。

2.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，步骤3.1中，渗透率、非Darcy流β因子、加速度系数的参考值k_r、β_r、c_ar、幂指数n_β和n_c的可能取值范围是：

k_r∈[k_r1,k_r2]，    （3-1）

β_r∈[β_r1,β_r2]，    （3-2）

c_ar∈[c_ar1,c_ar2]，    （3-3）

n_β∈[n_β1,n_β2]，    （3-4）

n_c∈[n_c1,n_c2]。    （3-5）

3.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，步骤3.3中所述遗传算法的个体基因型，相应的表现型为：

$k_r^j=k_{r1}{\left[\exp \frac{\ln \frac{k_{r2}}{k_{r1}}}{2^{l_1}-1}\right]}^j,j=\underset{i=1}{\overset{l_1}{\mathrm{\Σ}}}{2}^i{I}_{1i},j=\mathrm{0,1},...,2^{l_1}-1---(3-6)$

${\mathrm{\β}}_r^j={\mathrm{\β}}_{r1}{\left[\exp \frac{\ln \frac{{\mathrm{\β}}_{r2}}{{\mathrm{\β}}_{r1}}}{2^{l_2}-1}\right]}^j,j=\underset{i=1}{\overset{l_2}{\mathrm{\Σ}}}{2}^i{I}_{2i},j=\mathrm{0,1},...,2^{l_2}-1---(3-7)$

$c_{\mathrm{ar}}^j=c_{\mathrm{ar}1}{\left[\exp \frac{\ln \frac{c_{\mathrm{ar}2}}{c_{\mathrm{ar}1}}}{2^{l_3}-1}\right]}^j,j=\underset{i=1}{\overset{l_3}{\mathrm{\Σ}}}{2}^i{I}_{3i},j=\mathrm{0,1},...,2^{l_3}-1---(3-8)$

$n_b^j=n_{b1}[1+\frac{j}{2^{l_4}-1}\frac{n_{b2}-n_{b1}}{n_{b1}}],j=\underset{i=1}{\overset{l_4}{\mathrm{\Σ}}}{2}^i{I}_{4i},j=\mathrm{0,1},...,2^{l_4}-1---(3-9)$

$n_c^j=n_{c1}[1+\frac{j}{2^{l_5}-1}\frac{n_{c2}-n_{c1}}{n_{c1}}],j=\underset{i=1}{\overset{l_5}{\mathrm{\Σ}}}{2}^i{I}_{5i},j=\mathrm{0,1},...,2^{l_5}-1---(3-10)$

4.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，步骤3.4.4中所述的个体基因的表现型为：

$k_r^j=k_{r1}\exp \left(\frac{j}{2^{l_1}-1}\ln \frac{k_{r2}}{k_{r1}}\right),j=\underset{i=1}{\overset{l_1}{\mathrm{\Σ}}}{2}^i{I}_{1i},j=\mathrm{0,1},...,2^{l_1}-1---(3-11)$

${\mathrm{\β}}_r^j={\mathrm{\β}}_{r1}\exp \left(\frac{j}{2^{l_2}-1}\ln \frac{{\mathrm{\β}}_{r2}}{{\mathrm{\β}}_{r1}}\right),j=\underset{i=1}{\overset{l_2}{\mathrm{\Σ}}}{2}^i{I}_{2i},j=\mathrm{0,1},...,2^{l_2}-1---(3-12)$

$c_{\mathrm{ar}}^j=c_{\mathrm{ar}1}\exp \frac{j}{2^{l_3}-1}\ln \frac{c_{\mathrm{ar}2}}{c_{\mathrm{ar}1}},j=\underset{i=1}{\overset{l_3}{\mathrm{\Σ}}}{2}^i{I}_{3i},j=\mathrm{0,1},...,2^{l_3}-1---(3-13)$

$n_b^j=n_{b1}[1+\frac{j}{2^{l_4}-1}\frac{n_{b2}-n_{b1}}{n_{b1}}],j=\underset{i=1}{\overset{l_4}{\mathrm{\Σ}}}{2}^i{I}_{4i},j=\mathrm{0,1},...,2^{l_4}-1---(3-14)$

$n_c^j=n_{c1}[1+\frac{j}{2^{l_5}-1}\frac{n_{c2}-n_{c1}}{n_{c1}}],j=\underset{i=1}{\overset{l_5}{\mathrm{\Σ}}}{2}^i{I}_{5i},j=\mathrm{0,1},...,2^{l_5}-1---(3-15)$

5.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，步骤4.2中所述的动量守恒方程的外部函数为：

$\frac{\∂V}{\∂t}=-\frac{1}{{\mathrm{\ρc}}_a}(G_p+\frac{\mathrm{\μ}}{k}V+\mathrm{\ρ\β}{V}^2)---(4-1)$

其中，ρ为水的质量密度，β为水的动力粘度；其中，t∈[t_k,t_k+2]时刻的渗透率、非Darcy流β因子、加速度系数和压力梯度按全区间不等距三节点Lagrange插值计算如下：

$k=k^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+k^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+k^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---(4-2)$

$\mathrm{\β}={\mathrm{\β}}^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+{\mathrm{\β}}^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+{\mathrm{\β}}^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---(4-3)$

$c_a=c_a^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+c_a^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+c_k^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}---(4-4)$

$G_p=G_p^{t_k}\frac{(t-t_{k+1})(t-t_{k+2})}{(t_k-t_{k+1})(t_k-t_{k+2})}+G_p^{t_{k+1}}\frac{(t-t_k)(t-t_{k+2})}{(t_{k+1}-t_k)(t_{k+1}-t_{k+2})}+G_p^{t_{k+2}}\frac{(t-t_k)(t-t_{k+1})}{(t_{k+2}-t_k)(t_{k+2}-t_{k+1})}.---(4-5)$

6.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，步骤4.4中所述的E_rr的计算公式为： $E_{\mathrm{rr}}=\frac{1}{N}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(1-\frac{V_i^{*}}{V_i})}^2.---(4-6)$

7.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，步骤4.5中所述的计算种群中个体的适应度的公式为：

$\mathrm{fitn}\left(i\right)=\frac{1}{{\mathrm{err}}^{1/4}},i=\mathrm{1,2},...,k_{\mathrm{group}}.---(4-7)$

8.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，步骤5.1中所述的具有交配权的k_cop个体的基因型为：

$P_{\mathrm{cop}}=\left\{\mathrm{chromosome}\left(i_1\right)\right|\mathrm{chromosome}\left(i_1\right)=I_1^{i_1}{I}_2^{i_1}...I_m^{i_1},i=\mathrm{1,2},...,k_{\mathrm{cop}}\}.---(5-1)$

9.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，经步骤5.2进行交叉运算后得到个体集合

经步骤5.3对中每一个体进行变异运算得到

新一代种群： $\mathrm{NEW}\_\mathrm{pop}=\{I_1^i{I}_2^i...I_m^i|i=\mathrm{1,2},...,k_{\mathrm{group}}\}.---(5-2)$

10.根据权利要求1中所述的一种变质量破碎岩石渗透性参量计算方法，其特征在于，步骤5.5中，如果新一代群体中个体适应度的最大值fitn_max大于等于预先设定的数值s，即fitn_max≥s则停止繁殖；如果fitn_max＜s，则继续进行选择、交叉、变异运算，直到fitn_max≥s得以满足或繁殖代数等于某一事先设定的值T，则停止繁殖。

Images