一類空氣溫度的隨機(jī)微分方程模型

徐明明， 田德生

(湖北工業(yè)大學(xué)理學(xué)院， 湖北 武漢 430068)

對(duì)于溫度問(wèn)題的研究，主要有氣象學(xué)和統(tǒng)計(jì)學(xué)兩類方法。氣象學(xué)的方法是通過(guò)對(duì)氣象云圖以及其他氣象資料進(jìn)行分析，綜合太陽(yáng)輻射、大氣環(huán)流等因素進(jìn)行分析預(yù)測(cè)。統(tǒng)計(jì)學(xué)的方法則通過(guò)收集歷史數(shù)據(jù)，然后從這些數(shù)據(jù)中找出規(guī)律。統(tǒng)計(jì)學(xué)研究的常用方法和模型有最優(yōu)子集回歸，均生函數(shù)模型[2]等。隨機(jī)微分方程自20世紀(jì)中葉發(fā)展起來(lái)到現(xiàn)在，其理論得到了快速發(fā)展，研究成果已經(jīng)在金融、物理、生物、控制論、信息學(xué)等學(xué)科都有著廣泛的應(yīng)用，如金融經(jīng)濟(jì)學(xué)中的短期利率模型和二叉模型[1]、物理學(xué)中的布朗粒子逃逸問(wèn)題等。本文擬研究全球氣溫的變化問(wèn)題，根據(jù)氣溫變化的特點(diǎn)和大趨勢(shì)，嘗試建立相應(yīng)的隨機(jī)微分方程模型，并對(duì)其進(jìn)行分析。

在1850年之前，地球的平均溫度是相對(duì)穩(wěn)定的。但在1860年之后，尤其是進(jìn)入了工業(yè)社會(huì)以后，人類大量使用煤、石油等礦物燃料，排放出大量的CO2等溫室氣體。溫室氣體的排放是導(dǎo)致全球氣候變暖的主要因素。從氣象學(xué)的角度來(lái)說(shuō)，溫度的變化主要受太陽(yáng)照射、地理位置、海拔高度等自然因素的影響，同時(shí)，氣溫也受到人類活動(dòng)因素的影響，如近幾十年溫室氣體的大量排放導(dǎo)致溫度的持續(xù)上升，之后在各國(guó)的共同努力下低碳減排，尤其是聯(lián)合國(guó)氣候變化框架公約和《巴黎協(xié)定》等政策性文件的簽訂，使得全球氣溫的增長(zhǎng)變慢。這說(shuō)明全球氣溫的變化還是有規(guī)律可循的。

全球氣溫的頻繁變化對(duì)生態(tài)環(huán)境和動(dòng)植物的生存，乃至人類的生存都有較大的影響。在歷史的長(zhǎng)河中，全球氣溫就是一個(gè)時(shí)間序列，其變化既受到確定因素的影響(表現(xiàn)為增減趨勢(shì))，又受到不確定因素的影響(表現(xiàn)為隨機(jī)的波動(dòng)性)。這種表現(xiàn)在理論層面上，與隨機(jī)微分方程解的表現(xiàn)有些相似。因此，本文通過(guò)研究全球氣溫變化特點(diǎn)，建立關(guān)于氣溫的隨機(jī)微分方程模型。

1 預(yù)備知識(shí)

1.1 維納過(guò)程和伊藤微積分

維納過(guò)程是一個(gè)無(wú)后效性的獨(dú)立增量過(guò)程，其在大于t時(shí)刻的概率特性只與t時(shí)刻所處的狀態(tài)有關(guān)，而與t時(shí)刻之前的狀態(tài)無(wú)關(guān)。維納過(guò)程也稱為布朗運(yùn)動(dòng)(Brown)過(guò)程。我們把滿足以下條件的隨機(jī)過(guò)程稱為維納過(guò)程或Brown運(yùn)動(dòng)：

1)B(0)=0；

2)B(t),t≥0是獨(dú)立增量過(guò)程且具有平穩(wěn)增量；

3)對(duì)于任意的t>0,B(t)是均值為0，方差為δ2t的正態(tài)隨機(jī)分量。

特別地，當(dāng)δ2=1時(shí)，我們稱這個(gè)過(guò)程為標(biāo)準(zhǔn)維納過(guò)程[3]或標(biāo)準(zhǔn)布朗運(yùn)動(dòng)。

引理1 (重對(duì)數(shù)律[4]) 對(duì)于維納過(guò)程B(t)有

其中，sup表示上確界，a.s.表示幾乎必然。

伊藤(Ito)微積分是由日本學(xué)者伊藤清在1951年首次提出的，其目的是通過(guò)引入伊藤微積分解決隨機(jī)過(guò)程微積分問(wèn)題。伊藤微分公式的計(jì)算規(guī)則如下：

dt·dt=dt·dB(t)=

dB(t)·dt=0,dB(t)·dB(t)=dt

考慮隨機(jī)微分方程

dX(t)=f(t,X(t))dt+g(t,X(t))dB(t)

(1)

其中：f(x,y),g(x,y)是二元連續(xù)函數(shù)，B(t)是一維維納過(guò)程。解的形式的隨機(jī)過(guò)程如下：

此時(shí)方程的解X(t)稱為Ito過(guò)程，或擴(kuò)散過(guò)程。

引理2 (Ito公式[5])設(shè)ξ為Ito過(guò)程，即

dξ=b(t)ξdt+σ(t)ξdB(t)

如果實(shí)值函數(shù)h(t,x)對(duì)x二次連續(xù)可微。此時(shí)，η(t)=h(t,ε)也是Ito過(guò)程。有

(2)

且根據(jù)伊藤微分計(jì)算規(guī)則有(dξ)2=σ2(t)dt。

假設(shè)函數(shù)H(t)是一個(gè)左連續(xù)(適應(yīng)的)局部的有界過(guò)程。令Δ>0,N=t/Δ,ti=iΔ,0≤i≤Δ，此時(shí)從0到時(shí)間t，函數(shù)H(t)關(guān)于維納過(guò)程B(t)的Ito積分是一個(gè)隨機(jī)變量

由隨機(jī)微分方程的相關(guān)知識(shí)可證，該極限依概率收斂。且該隨機(jī)變量的期望滿足

1.2 平凡解的穩(wěn)定性

考慮自治的隨機(jī)微分方程問(wèn)題[6]

(3)

其中b(x)，σ(x)是兩個(gè)連續(xù)可微的函數(shù)，且滿足線性增長(zhǎng)條件。若同時(shí)滿足條件|b(x)|+|σ(x)|

1.3 隨機(jī)微分方程的數(shù)值計(jì)算

隨機(jī)微分方程的數(shù)值計(jì)算方法有很多，其中比較經(jīng)典的兩種分別是Euler法和Milstein法。這里主要介紹一下最簡(jiǎn)單也是運(yùn)用最廣泛的Euler法。

對(duì)方程(1)的解在區(qū)間[t0,T]上進(jìn)行離散化處理，然后利用Euler法構(gòu)造方程(1)的連續(xù)解，其Euler顯示迭代公式如下

X(tn+1)=X(tn)+f(tn,X(tn))(tn-tn-1)+g(tn,X(tn))(B(tn)-B(tn-1))

(4)

1.4 矩法估計(jì)

設(shè)總體ξ具有分布函數(shù)F(x,θ1,θ2,…,θk)，其中θ1,θ2,…,θk為未知參數(shù)，x1,x2,…,xn是取自總體ξ的一個(gè)樣本，又假設(shè)該總體的k階原點(diǎn)矩μk=Eξk存在，則有矩方程組

(5)

2 模型的建立與分析

2.1 氣溫變化的隨機(jī)微分方程模型

世界各地的人類活動(dòng)在不同時(shí)間段是不確定性的，導(dǎo)致全球氣溫改變的隨機(jī)性。一般地，當(dāng)前的氣溫值往往會(huì)影響下一時(shí)刻的溫度變化。以T(t)表示t時(shí)刻的全球大氣溫度，這是一個(gè)時(shí)間序列。全球氣溫的變化只與當(dāng)前時(shí)刻有關(guān)而與之前的時(shí)刻無(wú)關(guān)，即假設(shè)T(t)是一個(gè)維納過(guò)程。因此，可用隨機(jī)微分方程來(lái)描述全球氣溫模型。設(shè)f(t,T)表示人類活動(dòng)量引起的溫度T隨時(shí)間t變化率，于是，可得如下的隨機(jī)微分方程模型

(6)

其中：B(t)表示一維標(biāo)準(zhǔn)維納過(guò)程；g(t,T)表示擴(kuò)散系數(shù)。

2.2 氣溫變化的隨機(jī)指數(shù)增長(zhǎng)模型

收集1909年-2017年的全球平均溫度數(shù)據(jù)，單位為作圖1。

數(shù)據(jù)來(lái)源:美國(guó)宇航局戈達(dá)德研究所地球政策研究所圖 1 1909年-2017年的全球平均溫度

由圖1可見，溫度變化總體上呈現(xiàn)不斷波動(dòng)和增長(zhǎng)勢(shì)態(tài)。因此，結(jié)合式(6)，以下式來(lái)描述全球平均溫度模型

(7)

式中：r表示內(nèi)稟增長(zhǎng)率,α表示噪聲系數(shù)。記1909年為t=0時(shí)刻，該年的平均溫度為T0,即

T(0)=T0

(8)

又記方程(7)滿足初始條件(8)的解為T(t;T0)。

證明 方程(7)的系數(shù)連續(xù)可微且滿足線性增長(zhǎng)條件，因此方程(7)滿足初始條件(8)的解T(t;T0)存在且唯一。

對(duì)于式(7)，由引理2Ito公式有

根據(jù)Ito微分計(jì)算規(guī)則，結(jié)合式(7)，可得

(9)

積分式(9)，并由初始條件(8)得

即

(10)

由式(10)可知，當(dāng)T0>0時(shí)，恒有T(t;T0)>0,?t>0。

證畢。

T(t;T0)的期望反映了全球平均溫度的變化情況，其方差往往反映了溫度的波動(dòng)。下面分別討論T(t;T0)的期望ET(t;T0)和方差VarT(t;T0)。

定理2 對(duì)于方程(7)，有

ET(t;T0)=T0exp(rt),VarT(t;T0)=

證明方程(7)的解可寫為如下的隨機(jī)過(guò)程

對(duì)該隨機(jī)過(guò)程兩邊取期望，并結(jié)合引理3，有

等式兩邊對(duì)t求導(dǎo)可得

因此，對(duì)方程進(jìn)一步求解可得

ET(t;T0)=T0exp(rt)

對(duì)T2(t;T0)由引理2Ito公式可得

dT2(t;T0)=(2r+α2)T2(t;T0)dt+2αT2(t;T0)dB(t)

其隨機(jī)過(guò)程為

等式兩邊取期望有

同理，等式兩邊對(duì)t求導(dǎo)，再進(jìn)一步求解可得

證畢。

考慮方程(7)的參數(shù)估計(jì)。令lnT(t;T0)=x(t)，記lnT0=x0。 則式(9)變?yōu)?/p>

(11)

對(duì)于Δt>0，以等距時(shí)間點(diǎn)列0,Δt,2Δt,…,nΔt。 離散化方程(11)可得

其中εi(i=1,2,…,n)是獨(dú)立且均服從標(biāo)準(zhǔn)正態(tài)分布。因此由矩方程組(5)可設(shè)

解這個(gè)方程組，并且注意到有

因此

所以，有

其中，Ti表示第i(i=0,1,2,…,n)年的平均溫度。

2.3 全球氣溫的隨機(jī)阻滯增長(zhǎng)模型

由定理2知，當(dāng)t→∞時(shí)，有ET(t;T0)→∞(t→∞)， 即隨著時(shí)間的推移，溫度會(huì)無(wú)限增長(zhǎng)。顯然這與實(shí)際情況不符，因?yàn)樽匀毁Y源畢竟有限，有限的環(huán)境條件會(huì)對(duì)溫度的增長(zhǎng)起到阻滯作用，全球溫度的上升會(huì)逐漸變慢并最終會(huì)在某個(gè)溫度附近波動(dòng)?；诖耍⒁粋€(gè)關(guān)于溫度的隨機(jī)阻滯增長(zhǎng)模型

(12)

其中Tm表示環(huán)境容納量，即自然資源，環(huán)境條件所能容納的最高溫度。根據(jù)引理2Ito公式有

再將式(12)代入，可得

(13)

仍記方程(12)滿足初始條件(8)的解為T(t;T0)，對(duì)式(13)積分有

(14)

證明當(dāng)初始值T0>0時(shí)，根據(jù)式(14)，顯然有T(t;T0)>0。再對(duì)式(12)進(jìn)行變形，有

(15)

則根據(jù)式(15)有T(t;T0)

證畢。

3 數(shù)值模擬

受到自然資源等各種客觀因素的限制，環(huán)境容納量的值是一定的，因此，借助R語(yǔ)言軟件中的deSolve包對(duì)Tm值進(jìn)行估計(jì)，得到的結(jié)果為Tm=16.6224。對(duì)于方程(12)其Euler顯示迭代公式為

T(tn)=T(tn-1)+rT(tn-1)·

αT(tn-1)[B(tn)-B(tn-1)]

由此可以模擬得到含隨機(jī)擾動(dòng)項(xiàng)和不含隨機(jī)擾動(dòng)項(xiàng)的模型擬合效果(圖2)。

圖 2 模型擬合效果

圖2中，直線表示不含隨機(jī)擾動(dòng)項(xiàng)的阻滯增長(zhǎng)模型，曲線表示模型(11)所擬合的效果。由于內(nèi)稟增長(zhǎng)率r值較小，因此，不含隨機(jī)擾動(dòng)項(xiàng)的模型增長(zhǎng)緩慢，在短期內(nèi)圖形呈直線增長(zhǎng)。顯然本文擬合的隨機(jī)阻滯增長(zhǎng)模型更加合理。

進(jìn)一步，對(duì)全球平均溫度進(jìn)行了一定量的預(yù)測(cè)(圖3)。

圖 3 全球平均溫度預(yù)測(cè)

由圖3可以看出，在未來(lái)的很長(zhǎng)一段時(shí)間內(nèi)，全球平均溫度還是會(huì)呈現(xiàn)緩慢、波動(dòng)的上漲趨勢(shì)。通過(guò)數(shù)據(jù)擬合得到的環(huán)境容納量為16.6224℃，也就是說(shuō)最終在上升了2℃左右溫度將趨于穩(wěn)定。溫度首次達(dá)到這一數(shù)值的時(shí)間t=2395，也就是再過(guò)375年左右全球平均溫度才會(huì)趨于穩(wěn)定，達(dá)到環(huán)境容納值。

全球平均氣溫上升2℃也是我們?nèi)祟惪梢猿惺艿臉O限。2015年11月30日在巴黎舉辦的聯(lián)合國(guó)氣候大會(huì)指出：如果到21世紀(jì)末，人類還無(wú)法將全球平均氣溫增幅控制在2℃以下，會(huì)導(dǎo)致海平面上升、森林火災(zāi)、水資源缺乏等一系列災(zāi)難性后果。這個(gè)結(jié)果也與2014年聯(lián)合國(guó)政府間氣候變化專門委員會(huì)(IPCC)在德國(guó)柏林發(fā)布的第五次評(píng)估報(bào)告的研究結(jié)果相同。這說(shuō)明本文擬合出的環(huán)境容納量的值是準(zhǔn)確有效的。