王 昕，程希明

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

0 引言

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

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

1 模型構(gòu)建

1.1 翁氏模型

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

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

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

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

1.2 麥克斯韋模型

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

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

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

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

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

2 模型性質(zhì)

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

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

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

2.2 盛衰階段

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

3 參數(shù)估計

3.1 列文伯格-麥夸特算法

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

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

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

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

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

3.2 一元多項式回歸法

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

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

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

3.3 二元線性回歸法

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

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

4 算例分析

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

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

4.2 計算結(jié)果比較

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

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

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

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

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

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

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

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

4.4 預(yù)測數(shù)據(jù)對比

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

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

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

5 結(jié)論

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

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

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