0 引言

由于石油资源的特殊性和有限性，准确预测未来的石油产量能够为国家石油资源的战略性规划提出可行性建议，因此，石油产量的准确预测问题一直受到国内外学者的关注。自20世纪开始，学者们运用数学模型对石油供给能力进行量化分析：如Hubbert^[1-2]提出的Hubbert曲线，描述了石油、天然气等不可再生能源的产量变化曲线和供给规律，改变了世界石油科技专家们对不可再生能源供给能力线性增长的认识，并进一步提升了对经济发展与石油、天然气等资源供给关系的预期；翁文波^[3]提出了泊松旋回（Poisson Cycle）和逻辑斯谛旋回（Logistic Cycle）预测模型，并于1991年把泊松旋回模型更名为生命旋回（Life Cycle）模型^[4]，指出该模型是收敛的，仅限用于有限体系，如石油、天然气等不可再生的矿产资源等，并将该模型运用到全球原油产量基值的宏观预测中，获得了很好的预测效果。在此基础上，国内外学者对翁氏模型的推导及参数估计作了大量扩展研究：陈元千等^[5-7]提出了求解翁氏模型的线性试差法，并通过无因次处理油田实际开发数据与典型曲线的最佳拟合，从而得到模型参数值；张虎俊等^[8]提出了一种只需要进行2次简单线性回归的参数求解方法；乐平等^[9]利用峰值时间，将翁氏模型中的三参数转化为两参数，从而将拟合过程变为简单的二元线性回归；何俊等^[10]应用产量与时间存在的线性关系，进行线性回归，通过比较相关系数，确定出产量递减类型；赵林等^[11]提出求解翁氏模型的二元回归方法；田淑芳等^[12]用翁氏旋回模型对已探明石油储量与时间的关系进行拟合，进而对辽河油田石油探明储量增长趋势进行预测；张旭等^[13]借助最优化理论，提出了应用Gauss-Newton法进行模型求解参数的思路，并对实际油田产量进行了拟合；黄全华等^[14]根据油田年产油量统计数据，建立了改进无偏灰色模型，考虑了未来不确定因素的扰动，提高了预测精度。近年来，一些学者在不同数学理论的基础上，提出了新的预测模型：岳世俊等^[15]提出了枚举平衡法进行油藏数值模拟的模型初始化，用于提高数值模型的精确度；潘有军等^[16]使用灰色理论和多元线性回归法对油井压裂后的初期产能进行预测，比传统的直接多元线性回归法预测精度更高；刘美佳等^[17]结合物质平衡关系建立了一种产量随时间递减关系的新模型，并提出了“分段直线法”相对渗透率曲线拟合模型及“分段预测模型法”预测递减阶段的产量。

本次研究基于概率统计中的麦克斯韦分布，建立只有2个参数的新模型，来弥补广义翁氏模型在求解时，三参数线性试差法人为干扰因素较多，且计算量较大的不足，并对这一模型使用3种不同的参数估计方法进行求解，通过实证数据检验，并对比计算结果，使模型结果更加合理可靠。

1 模型构建 1.1 翁氏模型

翁文波^[3]提出，对于生命总量有限的许多体系，如非再生资源，可假设其产量在随时间t的变化过程中，正比于t的n次方函数，又随着t的负指数函数衰减，这一过程可以表示为

$ {Q_{\rm{t}}} = A{t^n}{{\rm{e}}^{ - t}},\;\;\;t \ge 0 $

(1)

式中：t表示油气田的开发时间，a；Q_t表示油气田产量，万t/a（油田）或亿m³/a（气田）；A表示模型参数，为正实数。式（1）为翁氏模型，但文献[3]中没有对该模型的建立给出推导过程。

陈元千等^[5-6]利用卡方分布和伽马分布，对翁氏模型进行了系统推导，得到如下形式的广义翁氏预测模型，并提出了求解模型的线性试差法

$ {Q_{\rm{t}}} = a{t^b}{{\rm{e}}^{ - \left( {\frac{t}{c}} \right)}} $

(2)

式（2）是一个含有3个参数的非线性模型，参数a，c为正实数，b为任意非负实数。在研究上述文献的基础上，发现线性试差法虽然易于理解，但预先设定不同的参数b值，用于得到参数a，c的值，这一步骤的人为因素较多，且不同的参数b值对a，c取值的影响较大。

1.2 麦克斯韦模型

作为连续型分布之一，麦克斯韦分布的概率密度函数可表示为

$ f(t) = \left\{ {\begin{array}{*{20}{l}} {\frac{4}{{b\sqrt {{\rm{ \mathsf{ π} }}b} }}{t^2}{{\rm{e}}^{ - \frac{{{t^2}}}{b}}}}&{t \ge 0}\\ 0&{t < 0} \end{array}} \right. $

(3)

式中：参数b > 0表示控制分布形态的常数。则有式（4）成立

$ \int\limits_0^\infty {f(t){\rm{d}}t} = 1 $

(4)

对于油田，投入开发时间t从0到∞的累计产量可视为可采储量N_R（万t/a），即

$ {N_{\rm{R}}} = \int\limits_0^\infty {Q(t){\rm{d}}t} $

(5)

式中：Q(t)表示第t年的石油产量，万t/a，若将式（5）转为预测油田产量的模型时，需要在式（5）右端乘以可采储量N_R。对于预测模型来说，N_R可被理解为由理论模型应用到实际问题时的模型转换系数。假设$\frac{Q(t)}{N_{\mathrm{R}}}=f(t)$，由式（5）可得

$ Q(t) = \frac{{4{N_{\rm{R}}}}}{{b\sqrt {{\rm{ \mathsf{ π} }}b} }}{t^2}{{\rm{e}}^{ - \frac{{{t^2}}}{b}}} $

(6)

令$A = \frac{{4{N_{\rm{R}}}}}{{b\sqrt {{\rm{ \mathit{ π} }}b} }}$，则式（6）可写为

$ Q(t) = A{t^2}{{\rm{e}}^{ - \frac{{{t^2}}}{b}}} $

(7)

其中A，b为待定参数，其值可以是任意的正实数，需由实际生产数据拟合而定。与广义翁氏模型相比，式（7）只有2个参数，避免了使用线性试差法求解时，需要预先设定1个参数值，才能求解另外2个参数值的人为影响。

2 模型性质 2.1 最高产量及其发生时间

由式（7）对生产时间t求导数，$\frac{{dQ(t)}}{{dt}} = 2Q(t) \cdot \left( {\frac{1}{t} - \frac{t}{b}} \right) $, 令$\frac{{dQ(t)}}{{dt}} = 0$可得最高年产量对应的生产时间

$ {t_{\rm{m}}} = \sqrt b $

(8)

将式（8）带入式（7）可得最高年产量（万t/a）和可采储量分别为

$ Q\left( {{t_{\rm{m}}}} \right) = At_{\rm{m}}^2{{\rm{e}}^{ - \frac{{{t^2}}}{b}}} $

(9)

$ {N_{\rm{R}}} = \frac{{Ab\sqrt {{\rm{ \mathsf{ π} }}b} }}{4} $

(10)

2.2 盛衰阶段

油气田的产量受到地下资源、世界经济、勘探开采技术等多个因素的影响，但产量大致的变化趋势可以分为4个阶段，分别是一般上升阶段、加速上升阶段、快速下降阶段和一般下降阶段，由式（7）可得

$ \frac{{{d^2}Q(t)}}{{d{t^2}}} = 2Q(t)\left( {\frac{1}{{{t^2}}} - \frac{5}{b} + \frac{{2{t^2}}}{{{b^2}}}} \right) $

(11)

$\frac{{{d^2}Q(t)}}{{d{t^2}}} = 0$有2个解，即

$ t = \frac{{\sqrt {(5 - \sqrt {17} )b} }}{2}\;或\;t = \frac{{\sqrt {(5 + \sqrt {17} )b} }}{2} $

(12)

由式（12）可知，当$0 < t \le \frac{{\sqrt {(5 - \sqrt {17} )b} }}{2}$时，产量处于一般上升阶段；当$\frac{{\sqrt {(5 - \sqrt {17} )b} }}{2} < t \le \sqrt b$时，产量处于加速上升阶段；当$\sqrt b < t \le \frac{{\sqrt {(5 + \sqrt {17} )b} }}{2}$时，产量处于快速下降阶段；当$\frac{{\sqrt {(5 + \sqrt {17} )b} }}{2} < t < + \infty $时，产量处于一般下降阶段。

3 参数估计 3.1 列文伯格-麦夸特算法

将式（7）视为非线性拟合问题，可采用最优化方法进行近似求解。将油田已经发生的历史产量数据视为向量Q ∈ Rⁿ，开发时间视为向量t ∈ R，其中Rⁿ表示n维向量空间。未知参数A、b构成向量$p = \left( {\begin{array}{*{20}{l}} A\\ b \end{array}} \right)$，则式（7）可表示为

$ \mathit{\boldsymbol{Q}} = f(\mathit{\boldsymbol{p}}) $

(13)

根据函数关系式f(·)与含噪声的历史实际产量数据Q，可采用列文伯格-麦夸特算法（LM，Levenberg-Marquardt）求解上述问题，得到参数A和b的数值近似解，具体步骤如下。

步骤1：取初始值p₀，置k = 0，λ₀ = 10^-3，γ = 10（也可以是其他大于1的数），计算终值控制常数${\varepsilon _{\rm{k}}} = \left\| {\mathit{\boldsymbol{Q}} - f\left( {{\mathit{\boldsymbol{p}}_{\rm{k}}}} \right)} \right\|$。

步骤2：计算Jacobi矩阵J_k，计算N_k = J_k^TJ_k + λ_k I，构造增量正规方程N_k · δ = J_k^T · ε_k，其中I为单位矩阵。

步骤3：求解增量正规方程得到向量δ_k

（1）如果$\left\| {\mathit{\boldsymbol{Q}} - f\left( {{\mathit{\boldsymbol{p}}_{\rm{k}}} + {\mathit{\boldsymbol{\delta }}_{\rm{k}}}} \right)} \right\| < {\varepsilon _{\rm{k}}}$，则令p_k+1= p_k + δ_k，若$\left\| {{\mathit{\boldsymbol{\delta }}_{\rm{k}}}} \right\| < {\varepsilon _{\rm{k}}}$，停止迭代，输出结果p_k+1；否则令${\lambda _{k + 1}} = \frac{{{\lambda _k}}}{\gamma }$，转到步骤2；

（2）如果$\left\| {\mathit{\boldsymbol{Q}} - f\left( {{\mathit{\boldsymbol{p}}_{\rm{k}}} + {\mathit{\boldsymbol{\delta }}_{\rm{k}}}} \right)} \right\| \ge {\varepsilon _{\rm{k}}}$，则令λ_k+1 = γλ_k，重新解增量正规方程得到向量δ_k，返回步骤1。

3.2 一元多项式回归法

借助于文献[6-7]中线性试差法的求解思想，将式（7）改写为

$ \frac{Q}{{{t^2}}} = A{{\rm{e}}^{ - \frac{{{t^2}}}{b}}} $

(14)

对式（14）两端取对数，可得

$ \ln \left( {\frac{Q}{{{t^2}}}} \right) = \ln A - \frac{{{t^2}}}{b} $

(15)

令$Y = \ln \left( {\frac{Q}{{{t^2}}}} \right), {X^2} = {t^2}$, ${\beta _0} = \ln A, {\beta _1} = - \frac{1}{b}$，可得一元二次多项式回归模型

$ Y = {\beta _0} + {\beta _1}{X^2} $

(16)

根据油田已经发生的历史产量数据及开发时间，求出二次曲线的参数β₀和β₁的值。

3.3 二元线性回归法

对式（7）两端取对数，可得

$ \ln Q = \ln A + 2\ln t - \frac{{{t^2}}}{b} $

(17)

令Y = ln Q，X₁ = ln t，X₂ = t²，β₀ = ln A，β₁ = 2，${\beta _2} = \frac{1}{b}$，得到如下二元线性回归模型

$ Y = {\beta _0} + {\beta _1}{X_1} + {\beta _2}{X_2} $

(18)

4 算例分析 4.1 数据选择

为了说明麦克斯韦模型的适用性，利用Norway石油公司的实际石油开发数据进行全面的历史数据拟合，并对Norway未来石油产量进行了合理预测。其中1971—2016年Norway石油实际产量历史数据来自《BP世界能源数据统计年鉴》。

4.2 计算结果比较

由表 1可看出，LM法拟合的残差平方和最大，说明该方法估计的参数值误差最大，但LM法拟合的最高产量生产年份更为准确；一元多项式回归与二元线性回归法所得参数值接近，二者所得均方差非常低，说明回归法拟合效果很好；决定系数表征回归方程在多大程度上解释了因变量的变化，或者说方程对观测值的拟合程度如何，从决定系数的角度来看，二元线性回归的拟合效果很好，虽然LM法的决定系数比一元多项式回归法的决定系数更高，但由于决定系数主要适用于线性回归，而LM法和一元多项式回归法都是非线性模型，因此其拟合效果不能用决定系数来评价。

4.3 拟合数据对比

根据上述3种方法，利用表 2中的实际产量数据，对1971—2016年Norway石油产量进行拟合，得到3种方法的拟合产量数据，并做出与实际产量拟合产量的对比图，结果见图 1和表 2。

由图 1可看出，Norway历史石油产量大致经历了初始一般上升阶段、加速上升阶段和快速下降阶段等阶段的变化。之所以出现这一变化，最重要的原因是油田储量是有限体系，单个油田的可开采量必定会被全部开采完毕；其次，由于石油开采的技术进步，石油开采的速度也必定会不断加快，越过最高产量之后，必定会出现快速下降。图 1中有2个波峰的原因较为复杂，其一可能是期间有新的油田被发现，其二可能是有新的开采技术投入，使得开采速度发生变化等等。实际历史累计产量和3种方法的拟合累计产量之间的对比如图 2所示。

4.4 预测数据对比

运用表 1中3种计算方法所得的参数，对未来Norway石油产量进行预测（图 3）。

在图 1-3中，2种回归方法的拟合产量曲线、累积产量曲线和预测产量曲线均基本重合。预测结果表明，自2002年起Norway的石油产量已进入下降阶段，如果Norway未来关于石油的进口战略不发生变化，且没有新的油田投入开发的话，至2075年Norway的石油产量将只剩下44万t，石油资源将面临枯竭。另外，随着Norway石油产量的下降，国家的石油净进口量一定会上升，且随着国家经济的发展，石油的消耗量还会进一步加大，该国家对原油净进口的依赖程度会越来越大。

5 结论

（1）使用麦克斯韦分布建立了石油产量预测模型，与广义翁氏模型不同的是，该模型中只有2个参数，估计方法简便，并应用非线性拟合、一元多项式回归拟合、二元线性回归拟合共3种方法进行模型的参数估计。

（2）在回归方法的计算中，可以不必选择搜集到的全部历史产量数据，而是从中选取与目标函数拟合关系比较好的数据段，可以得到更准确的参数值和拟合方程。

（3）3种参数估计方法的拟合效果与数据段的选取策略有关，方法并不存在根本性差异，各有优劣，在实践应用中可互为补充。