有限个平行一级反应成烃动力学模型参数的标定方法

摘要

有限个平行一级反应成烃动力学模型参数的标定方法，应用于油气资源评价技术领域。建立在成烃动力学原理基础之上用以描述有机质成烃过程，该标定方法能有效的将干酪根成油、干酪根成气以及油成气这几种差别较大，并且能对不同类型有机质的动力学过程区分开来讨论。另外，在模型标定的数学方法上进行了优化近似求解，一是摒弃了无限个平行一级反应所带来的标定精度提升有限但相对沉余的计算量，二是通过对动力学参数的近似求解实现了实验条件下无精确解的求取。有效地解决了对有机质成烃的复杂过程的描述，实现了对有机质生烃特征的定量、动态地描述。

主权项

1.有限个平行一级反应成烃动力学模型参数的标定方法，其特征在于：首先，建立成烃的化学动力学模型，成烃的化学动力学模型包括干酪根成油的化学动力学模型、干酪根成气的化学动力学模型和油裂解成气的化学动力学模型；然后，将观测参数或实验参数输入到建立好的成烃的化学动力学模型，对成烃的化学动力学模型进行标定；其次，采用构造目标函数和惩罚函数优化求取近似精确解，实现对有机质成烃过程的动力学描述；最后，将标定后的成烃的化学动力学模型，应用到探区的有机质成烃史的定量、动态评价；有限个平行一级反应成烃动力学模型，包括干酪根成油的化学动力学模型、干酪根成气的化学动力学模型、油裂解成气的化学动力学模型；A、干酪根成油的化学动力学模型设干酪根KEO成油过程由一系列NO个平行一级反应组成，每个反应对应的活化能为EO_i，单位：千焦/摩尔，指前因子AO_i，单位：每分钟，并设对应每一个反应的干酪根的原始潜量为XO_i0，KO_i为第i个干酪根成油反应的反应速率常数，每一个反应生成的油为O_i，每一个反应的生油量为XO_i，i=1，2，…，NO，即至时间t分钟，时$\frac{d{\mathrm{XO}}_i}{\mathrm{dt}}={\mathrm{KO}}_i({\mathrm{XO}}_{i0}-{\mathrm{XO}}_i)---\left(2\right)$由阿伦尼乌斯方程得到：${\mathrm{KO}}_i={\mathrm{AO}}_i\exp \left(\frac{-E{O}_i}{\mathrm{RT}}\right)---\left(3\right)$i=1，2，…，NO其中R为气体常数，R的值为8.31447kJ/mol·K；T为热力学温度，单位：开尔文；实验采用恒速升温，升温速率为D，单位：摄氏度/分钟；$\frac{\mathrm{dT}}{\mathrm{dt}}=D,$即$\mathrm{dt}=\frac{\mathrm{dT}}{D}---\left(4\right)$由(2)～(4)式得$\frac{d{\mathrm{XO}}_i}{{\mathrm{XO}}_{i0}-{\mathrm{XO}}_i}=\frac{{\mathrm{AO}}_i}{D}·\exp (-\frac{{\mathrm{EO}}_i}{\mathrm{RT}})\mathrm{dT}$将上式从T0→T积分，并注意到XO_i(T0)=0，XO_i(T)=XO_i得${\mathrm{XO}}_i={\mathrm{XO}}_{i0}(1-\exp (-{\int }_{T_0}^T\frac{{\mathrm{AO}}_i}{D}·\exp (-\frac{{\mathrm{EO}}_i}{\mathrm{RT}})\mathrm{dT}))---\left(5\right)$NO个平行反应的总生油量XO的相对百分比则为$\mathrm{XO}=\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_i=\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}({\mathrm{XO}}_{i0}(1-\exp (-{\int }_{T_0}^T\frac{{\mathrm{AO}}_i}{D}·\exp (-\frac{\mathrm{EOi}}{\mathrm{RT}})\mathrm{DT}))---\left(6\right)$B、干酪根成气的化学动力学模型同理，若设干酪根直接成气的反应由NG个平行反应组成，每个平行反应的活化能为EG_i，单位：千焦/摩尔，指前因子AG_i，单位：每分钟，初始潜量为XG_i0，单位：无量纲，D为升温速率，单位：摄氏度/分钟，R为气体常数，得随温度变化的T时期干酪根KEO直接生气量百分比数的计算公式为$\mathrm{XG}=\underset{i=1}{\overset{\mathrm{NG}}{\Sigma }}{\mathrm{XG}}_i=\underset{i=1}{\overset{N}{\Sigma }}({\mathrm{XG}}_{i0}(1-\exp (-{\int }_{T_0}^T\frac{{\mathrm{AG}}_i}{D}·\exp \left(1\frac{{\mathrm{EG}}_i}{D}\right)\mathrm{dT}))---\left(7\right)$与(6)式相比，(7)式仅仅是有关变量的副标不同；O表示油，G表示气；C、油裂解成气的化学动力学模型设油裂解成气的过程由NOG个平行反应组成，每一反应的活化能为EOG_i，单位：千焦/摩尔，指前因子为AOG_i，单位：每分钟，对应的原始潜量为XOG_i0，单位：无量纲，第i个干酪根成油反应的反应速率常数为KOG_i0，当反应进行至时间t，单位：分钟时，产气率，用占总反应量的百分数表示为XOG_i，单位：无量纲，则有$\frac{d{\mathrm{XOG}}_i}{\mathrm{dt}}={\mathrm{KOG}}_{i0}·({\mathrm{XOG}}_{i0}-{\mathrm{XOG}}_i)---\left(8\right)$其中${\mathrm{KOG}}_i={\mathrm{AOG}}_i·\exp (-\frac{{\mathrm{EOG}}_i}{\mathrm{RT}})---\left(9\right)$为第i个油裂解成气反应的反应速率常数；(9)代入(8)式整理后对时间积分，并注意到XOG_i(t=0)=0，XOG_i(t)=XOG_i，得到NOG个平行反应的总生气量为$\mathrm{XOG}=\underset{i=1}{\overset{\mathrm{NOG}}{\Sigma }}{\mathrm{XOG}}_i=\underset{i=1}{\overset{\mathrm{NOG}}{\Sigma }}{\mathrm{XOG}}_{i0}(1-\exp ({\int }_{T0}^T-({\mathrm{AOG}}_i·\exp (-\frac{{\mathrm{EOG}}_i}{\mathrm{RT}})\mathrm{dT}))---\left(10\right)$同样，式中的R为气体常数，T为热力学温度，开尔文；最后，可以通过成烃的化学动力学模型的标定，实现应用探区的有机质成烃史的定量、动态评价；对成烃的化学动力学模型进行标定的步骤是：以干酪根成油的化学动力学模型的标定进行说明，干酪根成气的化学动力学模型和油裂解成气的化学动力学模型标定的方法相同；步骤1，构造目标函数设在某一升温速率l，达到某一温度j时由实验所测得的产油率为XO1_lj，单位：百分比，在相同的条件下，如果存在某一组EO_i，EO_i表示第i个活化能的大小，单位：千焦/摩尔，AO_i，表示第i个活化能对应的指前因子的大小，单位：每分钟，XO_i0，第i个反应的干酪根的原始潜量，单位：无量纲，XO_i0的取值使对所有的l，j都有XO_lj-XO_lj=0，则该组EO_i，AO_i，XO_i0，即为所求；但由于实验误差方面的原因，这实际上是不可能的；因此，只能求使XO_lj-XO_lj尽量小的EO_i，AO_i、XO_i0的取值；为此，构造目标函数是：$Q({\mathrm{EO}}_i,{\mathrm{AO}}_i,{\mathrm{XO}}_{i0})=\underset{l=1}{\overset{L0}{\Sigma }}\underset{j=1}{\overset{J0}{\Sigma }}{\left(\frac{{\mathrm{XO}1}_{\mathrm{lj}}-{\mathrm{XO}}_{\mathrm{lj}}}{{\mathrm{XP}1}_{\mathrm{lj}}}\right)}^2---\left(11\right)$式中L0为不同升温速率实验的数目，J0为从一条实验曲线上的采样点数；平行反应的数目取值越大，就越能包括干酪根成油的所有反应类型，因而应该越精确，但由于此时模型标定及随后模型应用时的计算量太大，难以实用化，而且，实际标定的过程中发现，虽然模型对实验条件下热模拟时间，单位：分钟，恒温温度，单位，摄氏度或定升温速率，单位，摄氏度每分钟，烃类产率，单位，百分比参数的拟合程度一般随平行反应的细分而改善，但平行反应的数目达到一定程度之后，拟合程度的改善已不明显；因此，只需用有限个具有一定间隔的平行反应；这样，由于EO_i通过确定平行反应的活化能的分布范围和相邻平行反应的活化能间隔而求解，则(11)式化为：$Q({\mathrm{AO}}_i,{\mathrm{XO}}_{i0})=\underset{l=1}{\overset{L0}{\Sigma }}\underset{j=1}{\overset{J0}{\Sigma }}{\left(\frac{{\mathrm{XO}1}_{\mathrm{lj}}-{\mathrm{XO}}_{\mathrm{lj}}}{{\mathrm{XO}1}_{\mathrm{lj}}}\right)}^2---\left(12\right)$另外，注意到干酪根成油模型$\mathrm{XO}=\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_i=\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}\left({\mathrm{XO}}_{i0}(1-\exp (-{\int }_{T_0}^T\frac{{\mathrm{AO}}_i}{D}·\exp (-\frac{\mathrm{EOi}}{\mathrm{RT}})\mathrm{DT})\right)---\left(6\right)$式中的AO_i、XO_i0，应满足：其中，(6)和(13)中NO表示平行反应活化能的个数；这样，模型(6)的确定，即动力学参数的求取，问题就化为求非负的目标函数(12)在满足约束条件(13)时的极小点问题；步骤2，构造惩罚函数上述含有约束条件的极小值的求解问题比较复杂，因为除了要使目标函数值逐渐下降之外，还要注意解的可行性，即看解是否处于约束条件所限定的范围之内，通过构建带放大系数的惩罚项，达到提升惩罚项约束力，使得惩罚项的收敛速度加快，提升标定的效率与精度，这里采用惩罚函数法将有约束极值问题(12)、(13)化为无约束极值问题，其步骤如下：A：对AO_i＞0，这一约束条件，化为无约束极值问题有即C₁(AO_i)＝[min(0,AO_i)]²放大系数$B_1=e^{-{\mathrm{AO}}_i}$$G_1=B_1·C_1=e^{-{\mathrm{AO}}_i}·{\left[\min (0,{\mathrm{AO}}_i)\right]}^2$B：对1≥XO_i0≥0，这一约束条件，化为无约束极值问题有放大系数即$G_2\left({\mathrm{XO}}_{\mathrm{io}}\right)=C_2\left({\mathrm{XO}}_{\mathrm{io}}\right)·B_2\left({\mathrm{XO}}_{\mathrm{io}}\right)=e^{-{\mathrm{xo}}_{\mathrm{io}}}·{\left[\min (0,{\mathrm{XO}}_{\mathrm{io}})\right]}^2+e^{{\mathrm{xo}}_{\mathrm{io}}-1}·{\left[\min (\mathrm{0,1}-{\mathrm{XO}}_{\mathrm{io}})\right]}^2$C：对这一约束条件，化为无约束极值问题有放大系数即$G_3\left({\mathrm{XO}}_{i0}\right)=C_3\left({\mathrm{XO}}_{\mathrm{io}}\right)·B_3\left({\mathrm{XO}}_{\mathrm{io}}\right)=e^{\frac{ϵ-|1-\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{\mathrm{io}}|}{ϵ}}·{\left[\min (0,\frac{ϵ-|1-\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{i0}|}{ϵ})\right]}^2$D：惩罚项的构建$G({\mathrm{XO}}_{i0},{\mathrm{AO}}_i)=G_1+G_2+G_3$$=e^{-{\mathrm{AO}}_i}·{\left[\min (0,{\mathrm{AO}}_i)\right]}^2+e^{-{\mathrm{xo}}_{\mathrm{io}}}·{\left[\min (0,{\mathrm{XO}}_{\mathrm{io}})\right]}^2$$+e^{{\mathrm{xo}}_{\mathrm{io}}-1}·{\left[\min (\mathrm{0,1}-{\mathrm{XO}}_{\mathrm{io}})\right]}^2+e^{\frac{ϵ-|1-\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{\mathrm{io}}|}{ϵ}}·{\left[\min (0,\frac{ϵ|1-\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{i0}|}{ϵ})\right]}^2---\left(14\right)$取一个充分大的正整数R1，由(12)和(14)式构造出惩罚函数F(AO_i,XO_i0)＝Q(AO_i,XO_i0)+R1·G(XO_i0,AO_i)        (15)如果所求出的极小点超出约束条件，则逐渐增大R1，当R1充分大时，(14)式的极小解即为目标函数(12)式的极小解，这样就将有约束极值问题化为相对容易求解的无约束极值问题；步骤3，求一阶偏导函数极小值存在的必要条件是：函数式(15)的一阶偏导数为0；A，先对目标函数求偏导：$\frac{\partial Q}{\partial {\mathrm{AO}}_m}=\underset{l=1}{\overset{L0}{\Sigma }}\underset{j=1}{\overset{J0}{\Sigma }}(-2\frac{{\mathrm{XO}1}_{\mathrm{lj}}-{\mathrm{XO}}_{\mathrm{lj}}}{{\mathrm{XO}1}_{\mathrm{lj}}^2}·\frac{\partial {\mathrm{XO}}_{\mathrm{lj}}}{\partial {\mathrm{AO}}_m})$其中：$\frac{\partial {\mathrm{XO}}_{\mathrm{lj}}}{\partial {\mathrm{AO}}_m}=\frac{\partial \underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}\left({\mathrm{XO}}_{i0}(1-\exp (-{\int }_{T_0}^T\frac{{\mathrm{QO}}_i}{D_l}\exp (-\frac{{\mathrm{EO}}_i}{\mathrm{RT}})\mathrm{dT}))\right)}{\partial {\mathrm{AO}}_m}$$\frac{\partial Q}{\partial {\mathrm{XO}}_{m0}}=\underset{;=1}{\overset{L0}{\Sigma }}\underset{j=1}{\overset{J0}{\Sigma }}(-2\frac{{\mathrm{XO}1}_{\mathrm{lj}}-{\mathrm{XO}}_{\mathrm{lj}}}{{\mathrm{XO}1}_{\mathrm{lj}}^2}·\frac{\partial {\mathrm{XO}}_{\mathrm{lj}}}{\partial {\mathrm{XO}}_{m0}})$$=\underset{l=1}{\overset{L0}{\Sigma }}\underset{j=1}{\overset{J0}{\Sigma }}(-2\frac{{\mathrm{XO}1}_{\mathrm{lj}}-{\mathrm{XO}}_{\mathrm{lj}}}{{\mathrm{XO}1}_{\mathrm{lj}}^2}·(1-\exp (-{\int }_{T_0}^T\frac{{\mathrm{AO}}_m}{D_l}\exp (-\frac{{\mathrm{EO}}_m}{\mathrm{RT}})\mathrm{dT})))$这里m=1，2，3，…，NO；B，对惩罚项求偏导：$\frac{\partial G}{\partial {\mathrm{AO}}_m}=2·e^{-{\mathrm{AO}}_m}·\min (0,{\mathrm{AO}}_m)-\min {(0,{\mathrm{AO}}_m)}^2·e^{-{\mathrm{AO}}_m}$$\frac{\partial G}{\partial {\mathrm{XO}}_{\mathrm{mo}}}=2·\min (0,{\mathrm{XO}}_{\mathrm{mo}})·e^{-{\mathrm{xo}}_{\mathrm{io}}}-e^{-{\mathrm{xo}}_{\mathrm{io}}}·{\left[\min (0,{\mathrm{XO}}_{\mathrm{io}})\right]}^2-2·e^{{\mathrm{xo}}_{\mathrm{io}}-1}·\min (\mathrm{0,1}-{\mathrm{XO}}_{\mathrm{mo}})$$+e^{{\mathrm{xo}}_{\mathrm{io}}-1}·{\left[\min (\mathrm{0,1}-x_{\mathrm{io}})\right]}^2+2·\min (0,\frac{ϵ-|1-\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{i0}}{ϵ})·\mathrm{FN}(\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{i0}-1)·e^{\frac{ϵ-|1-\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{i0}|}{ϵ}}$$+e^{\frac{ϵ-|1-\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{i0}|}{ϵ}}·\min {(0,\frac{ϵ-|1-\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{i0}}{ϵ})}^2·\mathrm{FN}(\underset{i=1}{\overset{\mathrm{NO}}{\Sigma }}{\mathrm{XO}}_{i0}-1)$这里FN为取括号里表达式的符号，即理论上讲，极小点处应有$\left\{\begin{array}{c}\frac{\partial F({\mathrm{AO}}_i,{\mathrm{XO}}_{i0})}{\partial {\mathrm{AO}}_m}=\frac{\partial Q}{\partial {\mathrm{AO}}_m}+R1·\frac{\partial G}{\partial {\mathrm{AO}}_m}=0\\ \frac{\partial F({\mathrm{AO}}_i,{\mathrm{XO}}_{i0})}{\partial {\mathrm{XO}}_{m0}}=\frac{\partial Q}{\partial {\mathrm{XO}}_{m0}}+R1·\frac{\partial G}{\partial {\mathrm{XO}}_{m0}}=0\\ m=\mathrm{1,2},···,\mathrm{NO}\end{array}\right.$因此，精确求出方程组(15)，即求解出2×NO个方程式，2×NO个待定变量，则得出若干可能的极小点，达到求解2×NO个待定动力学参数AO_i、XO_i0，标定模型${\mathrm{XO}}_i={\mathrm{XO}}_{i0}(1-\exp (-{\int }_{T_0}^T\frac{{\mathrm{AO}}_i}{D}·\exp (-\frac{{\mathrm{EO}}_i}{\mathrm{RT}})\mathrm{dT}))---\left(5\right)$的目的；但是对方程组(15)这样复杂的非多项式函数，不可能求出其精确解，但能求出其近似解；步骤4，近似极小点的求取对无约束值问题(15)的求解，采用收敛速度较快而又无需计算烦琐的二阶导数矩阵及其逆矩阵的变尺度法来进行优化计算。