王 昕，程希明

（北京信息科技大學(xué)理學(xué)院，北京 100192）

0 引言

由于石油資源的特殊性和有限性，準(zhǔn)確預(yù)測(cè)未來(lái)的石油產(chǎn)量能夠?yàn)閲?guó)家石油資源的戰(zhàn)略性規(guī)劃提出可行性建議，因此，石油產(chǎn)量的準(zhǔn)確預(yù)測(cè)問(wèn)題一直受到國(guó)內(nèi)外學(xué)者的關(guān)注。自20 世紀(jì)開(kāi)始，學(xué)者們運(yùn)用數(shù)學(xué)模型對(duì)石油供給能力進(jìn)行量化分析:如Hubbert［1-2］提出的Hubbert 曲線，描述了石油、天然氣等不可再生能源的產(chǎn)量變化曲線和供給規(guī)律，改變了世界石油科技專家們對(duì)不可再生能源供給能力線性增長(zhǎng)的認(rèn)識(shí)，并進(jìn)一步提升了對(duì)經(jīng)濟(jì)發(fā)展與石油、天然氣等資源供給關(guān)系的預(yù)期；翁文波［3］提出了泊松旋回（Poisson Cycle）和邏輯斯諦旋回（Logistic Cycle）預(yù)測(cè)模型，并于1991 年把泊松旋回模型更名為生命旋回（Life Cycle）模型［4］，指出該模型是收斂的，僅限用于有限體系，如石油、天然氣等不可再生的礦產(chǎn)資源等，并將該模型運(yùn)用到全球原油產(chǎn)量基值的宏觀預(yù)測(cè)中，獲得了很好的預(yù)測(cè)效果。在此基礎(chǔ)上，國(guó)內(nèi)外學(xué)者對(duì)翁氏模型的推導(dǎo)及參數(shù)估計(jì)作了大量擴(kuò)展研究:陳元千等［5-7］提出了求解翁氏模型的線性試差法，并通過(guò)無(wú)因次處理油田實(shí)際開(kāi)發(fā)數(shù)據(jù)與典型曲線的最佳擬合，從而得到模型參數(shù)值；張虎俊等［8］提出了一種只需要進(jìn)行2 次簡(jiǎn)單線性回歸的參數(shù)求解方法；樂(lè)平等［9］利用峰值時(shí)間，將翁氏模型中的三參數(shù)轉(zhuǎn)化為兩參數(shù)，從而將擬合過(guò)程變?yōu)楹?jiǎn)單的二元線性回歸；何俊等［10］應(yīng)用產(chǎn)量與時(shí)間存在的線性關(guān)系，進(jìn)行線性回歸，通過(guò)比較相關(guān)系數(shù)，確定出產(chǎn)量遞減類型；趙林等［11］提出求解翁氏模型的二元回歸方法；田淑芳等［12］用翁氏旋回模型對(duì)已探明石油儲(chǔ)量與時(shí)間的關(guān)系進(jìn)行擬合，進(jìn)而對(duì)遼河油田石油探明儲(chǔ)量增長(zhǎng)趨勢(shì)進(jìn)行預(yù)測(cè)；張旭等［13］借助最優(yōu)化理論，提出了應(yīng)用Gauss-Newton法進(jìn)行模型求解參數(shù)的思路，并對(duì)實(shí)際油田產(chǎn)量進(jìn)行了擬合；黃全華等［14］根據(jù)油田年產(chǎn)油量統(tǒng)計(jì)數(shù)據(jù)，建立了改進(jìn)無(wú)偏灰色模型，考慮了未來(lái)不確定因素的擾動(dòng)，提高了預(yù)測(cè)精度。近年來(lái)，一些學(xué)者在不同數(shù)學(xué)理論的基礎(chǔ)上，提出了新的預(yù)測(cè)模型:岳世俊等［15］提出了枚舉平衡法進(jìn)行油藏?cái)?shù)值模擬的模型初始化，用于提高數(shù)值模型的精確度；潘有軍等［16］使用灰色理論和多元線性回歸法對(duì)油井壓裂后的初期產(chǎn)能進(jìn)行預(yù)測(cè)，比傳統(tǒng)的直接多元線性回歸法預(yù)測(cè)精度更高；劉美佳等［17］結(jié)合物質(zhì)平衡關(guān)系建立了一種產(chǎn)量隨時(shí)間遞減關(guān)系的新模型，并提出了“分段直線法”相對(duì)滲透率曲線擬合模型及“分段預(yù)測(cè)模型法”預(yù)測(cè)遞減階段的產(chǎn)量。

本次研究基于概率統(tǒng)計(jì)中的麥克斯韋分布，建立只有2 個(gè)參數(shù)的新模型，來(lái)彌補(bǔ)廣義翁氏模型在求解時(shí)，三參數(shù)線性試差法人為干擾因素較多，且計(jì)算量較大的不足，并對(duì)這一模型使用3 種不同的參數(shù)估計(jì)方法進(jìn)行求解，通過(guò)實(shí)證數(shù)據(jù)檢驗(yàn)，并對(duì)比計(jì)算結(jié)果，使模型結(jié)果更加合理可靠。

1 模型構(gòu)建

1.1 翁氏模型

翁文波［3］提出，對(duì)于生命總量有限的許多體系，如非再生資源，可假設(shè)其產(chǎn)量在隨時(shí)間t的變化過(guò)程中，正比于t的n次方函數(shù)，又隨著t的負(fù)指數(shù)函數(shù)衰減，這一過(guò)程可以表示為

式中:t表示油氣田的開(kāi)發(fā)時(shí)間，a；Qt表示油氣田產(chǎn)量，萬(wàn)t/a（油田）或億m3/a（氣田）；A表示模型參數(shù)，為正實(shí)數(shù)。式（1）為翁氏模型，但文獻(xiàn)［3］中沒(méi)有對(duì)該模型的建立給出推導(dǎo)過(guò)程。

陳元千等［5-6］利用卡方分布和伽馬分布，對(duì)翁氏模型進(jìn)行了系統(tǒng)推導(dǎo)，得到如下形式的廣義翁氏預(yù)測(cè)模型，并提出了求解模型的線性試差法

式（2）是一個(gè)含有3 個(gè)參數(shù)的非線性模型，參數(shù)a，c為正實(shí)數(shù)，b為任意非負(fù)實(shí)數(shù)。在研究上述文獻(xiàn)的基礎(chǔ)上，發(fā)現(xiàn)線性試差法雖然易于理解，但預(yù)先設(shè)定不同的參數(shù)b值，用于得到參數(shù)a，c的值，這一步驟的人為因素較多，且不同的參數(shù)b值對(duì)a，c取值的影響較大。

1.2 麥克斯韋模型

作為連續(xù)型分布之一，麥克斯韋分布的概率密度函數(shù)可表示為

式中:參數(shù)b＞0 表示控制分布形態(tài)的常數(shù)。則有式（4）成立

對(duì)于油田，投入開(kāi)發(fā)時(shí)間t從0 到∞的累計(jì)產(chǎn)量可視為可采儲(chǔ)量NR（萬(wàn)t/a），即

式中:Q(t)表示第t年的石油產(chǎn)量，萬(wàn)t/a，若將式（5）轉(zhuǎn)為預(yù)測(cè)油田產(chǎn)量的模型時(shí)，需要在式（5）右端乘以可采儲(chǔ)量NR。對(duì)于預(yù)測(cè)模型來(lái)說(shuō)，NR可被理解為由理論模型應(yīng)用到實(shí)際問(wèn)題時(shí)的模型轉(zhuǎn)換系數(shù)。假設(shè)由式（5）可得

其中A，b為待定參數(shù)，其值可以是任意的正實(shí)數(shù)，需由實(shí)際生產(chǎn)數(shù)據(jù)擬合而定。與廣義翁氏模型相比，式（7）只有2 個(gè)參數(shù)，避免了使用線性試差法求解時(shí)，需要預(yù)先設(shè)定1 個(gè)參數(shù)值，才能求解另外2個(gè)參數(shù)值的人為影響。

2 模型性質(zhì)

2.1 最高產(chǎn)量及其發(fā)生時(shí)間

由式（7）對(duì)生產(chǎn)時(shí)間t求導(dǎo)數(shù)，可得最高年產(chǎn)量對(duì)應(yīng)的生產(chǎn)時(shí)間

將式（8）帶入式（7）可得最高年產(chǎn)量（萬(wàn)t/a）和可采儲(chǔ)量分別為

2.2 盛衰階段

油氣田的產(chǎn)量受到地下資源、世界經(jīng)濟(jì)、勘探開(kāi)采技術(shù)等多個(gè)因素的影響，但產(chǎn)量大致的變化趨勢(shì)可以分為4 個(gè)階段，分別是一般上升階段、加速上升階段、快速下降階段和一般下降階段，由式（7）可得

3 參數(shù)估計(jì)

3.1 列文伯格-麥夸特算法

將式（7）視為非線性擬合問(wèn)題，可采用最優(yōu)化方法進(jìn)行近似求解。將油田已經(jīng)發(fā)生的歷史產(chǎn)量數(shù)據(jù)視為向量Q∈Rn，開(kāi)發(fā)時(shí)間視為向量t∈R，其中Rn表示n維向量空間。未知參數(shù)A、b構(gòu)成向量，則式（7）可表示為

根據(jù)函數(shù)關(guān)系式f(?)與含噪聲的歷史實(shí)際產(chǎn)量數(shù)據(jù)Q，可采用列文伯格-麥夸特算法（LM，Levenberg-Marquardt）求解上述問(wèn)題，得到參數(shù)A和b的數(shù)值近似解，具體步驟如下。

步驟1:取初始值p0，置k=0，λ0=10-3，γ=10（也可以是其他大于1 的數(shù)），計(jì)算終值控制常數(shù)

步驟2:計(jì)算Jacobi 矩陣Jk，計(jì)算λkI，構(gòu)造增量正規(guī)方程，其中I為單位矩陣。

步驟3:求解增量正規(guī)方程得到向量δk

3.2 一元多項(xiàng)式回歸法

借助于文獻(xiàn)［6-7］中線性試差法的求解思想，將式（7）改寫(xiě)為

對(duì)式（14）兩端取對(duì)數(shù)，可得

根據(jù)油田已經(jīng)發(fā)生的歷史產(chǎn)量數(shù)據(jù)及開(kāi)發(fā)時(shí)間，求出二次曲線的參數(shù)β0和β1的值。

3.3 二元線性回歸法

對(duì)式（7）兩端取對(duì)數(shù)，可得

令Y=lnQ，X1=lnt，X2=t2，β0=lnA，β1=2，，得到如下二元線性回歸模型

4 算例分析

4.1 數(shù)據(jù)選擇

為了說(shuō)明麥克斯韋模型的適用性，利用Norway石油公司的實(shí)際石油開(kāi)發(fā)數(shù)據(jù)進(jìn)行全面的歷史數(shù)據(jù)擬合，并對(duì)Norway 未來(lái)石油產(chǎn)量進(jìn)行了合理預(yù)測(cè)。其中1971—2016 年Norway 石油實(shí)際產(chǎn)量歷史數(shù)據(jù)來(lái)自《BP 世界能源數(shù)據(jù)統(tǒng)計(jì)年鑒》。

4.2 計(jì)算結(jié)果比較

由表1 可看出，LM 法擬合的殘差平方和最大，說(shuō)明該方法估計(jì)的參數(shù)值誤差最大，但LM 法擬合的最高產(chǎn)量生產(chǎn)年份更為準(zhǔn)確；一元多項(xiàng)式回歸與二元線性回歸法所得參數(shù)值接近，二者所得均方差非常低，說(shuō)明回歸法擬合效果很好；決定系數(shù)表征回歸方程在多大程度上解釋了因變量的變化，或者說(shuō)方程對(duì)觀測(cè)值的擬合程度如何，從決定系數(shù)的角度來(lái)看，二元線性回歸的擬合效果很好，雖然LM法的決定系數(shù)比一元多項(xiàng)式回歸法的決定系數(shù)更高，但由于決定系數(shù)主要適用于線性回歸，而LM法和一元多項(xiàng)式回歸法都是非線性模型，因此其擬合效果不能用決定系數(shù)來(lái)評(píng)價(jià)。

表1 3 種參數(shù)方法估計(jì)結(jié)果Table 1 Estimation results of three methods

4.3 擬合數(shù)據(jù)對(duì)比

根據(jù)上述3 種方法，利用表2 中的實(shí)際產(chǎn)量數(shù)據(jù)，對(duì)1971—2016 年Norway 石油產(chǎn)量進(jìn)行擬合，得到3 種方法的擬合產(chǎn)量數(shù)據(jù)，并做出與實(shí)際產(chǎn)量擬合產(chǎn)量的對(duì)比圖，結(jié)果見(jiàn)圖1 和表2。

由圖1 可看出，Norway 歷史石油產(chǎn)量大致經(jīng)歷了初始一般上升階段、加速上升階段和快速下降階段等階段的變化。之所以出現(xiàn)這一變化，最重要的原因是油田儲(chǔ)量是有限體系，單個(gè)油田的可開(kāi)采量必定會(huì)被全部開(kāi)采完畢；其次，由于石油開(kāi)采的技術(shù)進(jìn)步，石油開(kāi)采的速度也必定會(huì)不斷加快，越過(guò)最高產(chǎn)量之后，必定會(huì)出現(xiàn)快速下降。圖1 中有2 個(gè)波峰的原因較為復(fù)雜，其一可能是期間有新的油田被發(fā)現(xiàn)，其二可能是有新的開(kāi)采技術(shù)投入，使得開(kāi)采速度發(fā)生變化等等。實(shí)際歷史累計(jì)產(chǎn)量和3種方法的擬合累計(jì)產(chǎn)量之間的對(duì)比如圖2 所示。

圖1 1971—2016 年Norway 石油實(shí)際產(chǎn)量與擬合產(chǎn)量Fig.1 Actual oil production and fitting oil production of Norway in 1971-2016

表2 1971—2016 年Norway 石油實(shí)際產(chǎn)量與擬合產(chǎn)量對(duì)比Table 2 Comparison between actual oil production and fitting oil production of Norway in 1971-2016

圖2 1971—2016 年Norway 石油累計(jì)實(shí)際產(chǎn)量和累計(jì)擬合產(chǎn)量Fig.2 Cumulative actual oil production and cumulative fitting oil production of Norway in 1971-2016

4.4 預(yù)測(cè)數(shù)據(jù)對(duì)比

運(yùn)用表1 中3 種計(jì)算方法所得的參數(shù)，對(duì)未來(lái)Norway 石油產(chǎn)量進(jìn)行預(yù)測(cè)（圖3）。

圖3 2017—2075 年Norway 石油產(chǎn)量預(yù)測(cè)Fig.3 Oil production forecasting of Norway in 2017-2075

在圖1—3 中，2 種回歸方法的擬合產(chǎn)量曲線、累積產(chǎn)量曲線和預(yù)測(cè)產(chǎn)量曲線均基本重合。預(yù)測(cè)結(jié)果表明，自2002 年起Norway 的石油產(chǎn)量已進(jìn)入下降階段，如果Norway 未來(lái)關(guān)于石油的進(jìn)口戰(zhàn)略不發(fā)生變化，且沒(méi)有新的油田投入開(kāi)發(fā)的話，至2075 年Norway 的石油產(chǎn)量將只剩下44 萬(wàn)t，石油資源將面臨枯竭。另外，隨著Norway 石油產(chǎn)量的下降，國(guó)家的石油凈進(jìn)口量一定會(huì)上升，且隨著國(guó)家經(jīng)濟(jì)的發(fā)展，石油的消耗量還會(huì)進(jìn)一步加大，該國(guó)家對(duì)原油凈進(jìn)口的依賴程度會(huì)越來(lái)越大。

5 結(jié)論

（1）使用麥克斯韋分布建立了石油產(chǎn)量預(yù)測(cè)模型，與廣義翁氏模型不同的是，該模型中只有2 個(gè)參數(shù)，估計(jì)方法簡(jiǎn)便，并應(yīng)用非線性擬合、一元多項(xiàng)式回歸擬合、二元線性回歸擬合共3 種方法進(jìn)行模型的參數(shù)估計(jì)。

（2）在回歸方法的計(jì)算中，可以不必選擇搜集到的全部歷史產(chǎn)量數(shù)據(jù)，而是從中選取與目標(biāo)函數(shù)擬合關(guān)系比較好的數(shù)據(jù)段，可以得到更準(zhǔn)確的參數(shù)值和擬合方程。

（3）3 種參數(shù)估計(jì)方法的擬合效果與數(shù)據(jù)段的選取策略有關(guān)，方法并不存在根本性差異，各有優(yōu)劣，在實(shí)踐應(yīng)用中可互為補(bǔ)充。