徐明明, 田德生
(湖北工業大學理學院, 湖北 武漢 430068)
對于溫度問題的研究,主要有氣象學和統計學兩類方法。氣象學的方法是通過對氣象云圖以及其他氣象資料進行分析,綜合太陽輻射、大氣環流等因素進行分析預測。統計學的方法則通過收集歷史數據,然后從這些數據中找出規律。統計學研究的常用方法和模型有最優子集回歸,均生函數模型[2]等。隨機微分方程自20世紀中葉發展起來到現在,其理論得到了快速發展,研究成果已經在金融、物理、生物、控制論、信息學等學科都有著廣泛的應用,如金融經濟學中的短期利率模型和二叉模型[1]、物理學中的布朗粒子逃逸問題等。本文擬研究全球氣溫的變化問題,根據氣溫變化的特點和大趨勢,嘗試建立相應的隨機微分方程模型,并對其進行分析。
在1850年之前,地球的平均溫度是相對穩定的。但在1860年之后,尤其是進入了工業社會以后,人類大量使用煤、石油等礦物燃料,排放出大量的CO2等溫室氣體。溫室氣體的排放是導致全球氣候變暖的主要因素。從氣象學的角度來說,溫度的變化主要受太陽照射、地理位置、海拔高度等自然因素的影響,同時,氣溫也受到人類活動因素的影響,如近幾十年溫室氣體的大量排放導致溫度的持續上升,之后在各國的共同努力下低碳減排,尤其是聯合國氣候變化框架公約和《巴黎協定》等政策性文件的簽訂,使得全球氣溫的增長變慢。這說明全球氣溫的變化還是有規律可循的。
全球氣溫的頻繁變化對生態環境和動植物的生存,乃至人類的生存都有較大的影響。在歷史的長河中,全球氣溫就是一個時間序列,其變化既受到確定因素的影響(表現為增減趨勢),又受到不確定因素的影響(表現為隨機的波動性)。這種表現在理論層面上,與隨機微分方程解的表現有些相似。因此,本文通過研究全球氣溫變化特點,建立關于氣溫的隨機微分方程模型。
維納過程是一個無后效性的獨立增量過程,其在大于t時刻的概率特性只與t時刻所處的狀態有關,而與t時刻之前的狀態無關。維納過程也稱為布朗運動(Brown)過程。我們把滿足以下條件的隨機過程稱為維納過程或Brown運動:
1)B(0)=0;
2)B(t),t≥0是獨立增量過程且具有平穩增量;
3)對于任意的t>0,B(t)是均值為0,方差為δ2t的正態隨機分量。
特別地,當δ2=1時,我們稱這個過程為標準維納過程[3]或標準布朗運動。
引理1 (重對數律[4]) 對于維納過程B(t)有
其中,sup表示上確界,a.s.表示幾乎必然。
伊藤(Ito)微積分是由日本學者伊藤清在1951年首次提出的,其目的是通過引入伊藤微積分解決隨機過程微積分問題。伊藤微分公式的計算規則如下:
dt·dt=dt·dB(t)=
dB(t)·dt=0,dB(t)·dB(t)=dt
考慮隨機微分方程
dX(t)=f(t,X(t))dt+g(t,X(t))dB(t)
(1)
其中:f(x,y),g(x,y)是二元連續函數,B(t)是一維維納過程。解的形式的隨機過程如下:
此時方程的解X(t)稱為Ito過程,或擴散過程。
引理2 (Ito公式[5])設ξ為Ito過程,即
dξ=b(t)ξdt+σ(t)ξdB(t)
如果實值函數h(t,x)對x二次連續可微。此時,η(t)=h(t,ε)也是Ito過程。有

(2)
且根據伊藤微分計算規則有(dξ)2=σ2(t)dt。
假設函數H(t)是一個左連續(適應的)局部的有界過程。令Δ>0,N=t/Δ,ti=iΔ,0≤i≤Δ,此時從0到時間t,函數H(t)關于維納過程B(t)的Ito積分是一個隨機變量
由隨機微分方程的相關知識可證,該極限依概率收斂。且該隨機變量的期望滿足
考慮自治的隨機微分方程問題[6]
(3)
其中b(x),σ(x)是兩個連續可微的函數,且滿足線性增長條件。若同時滿足條件|b(x)|+|σ(x)| 隨機微分方程的數值計算方法有很多,其中比較經典的兩種分別是Euler法和Milstein法。這里主要介紹一下最簡單也是運用最廣泛的Euler法。 對方程(1)的解在區間[t0,T]上進行離散化處理,然后利用Euler法構造方程(1)的連續解,其Euler顯示迭代公式如下 X(tn+1)=X(tn)+f(tn,X(tn))(tn-tn-1)+g(tn,X(tn))(B(tn)-B(tn-1)) (4) 設總體ξ具有分布函數F(x,θ1,θ2,…,θk),其中θ1,θ2,…,θk為未知參數,x1,x2,…,xn是取自總體ξ的一個樣本,又假設該總體的k階原點矩μk=Eξk存在,則有矩方程組 (5) 世界各地的人類活動在不同時間段是不確定性的,導致全球氣溫改變的隨機性。一般地,當前的氣溫值往往會影響下一時刻的溫度變化。以T(t)表示t時刻的全球大氣溫度,這是一個時間序列。全球氣溫的變化只與當前時刻有關而與之前的時刻無關,即假設T(t)是一個維納過程。因此,可用隨機微分方程來描述全球氣溫模型。設f(t,T)表示人類活動量引起的溫度T隨時間t變化率,于是,可得如下的隨機微分方程模型 (6) 其中:B(t)表示一維標準維納過程;g(t,T)表示擴散系數。 收集1909年-2017年的全球平均溫度數據,單位為作圖1。 數據來源:美國宇航局戈達德研究所地球政策研究所圖 1 1909年-2017年的全球平均溫度 由圖1可見,溫度變化總體上呈現不斷波動和增長勢態。因此,結合式(6),以下式來描述全球平均溫度模型 (7) 式中:r表示內稟增長率,α表示噪聲系數。記1909年為t=0時刻,該年的平均溫度為T0,即 T(0)=T0 (8) 又記方程(7)滿足初始條件(8)的解為T(t;T0)。 證明 方程(7)的系數連續可微且滿足線性增長條件,因此方程(7)滿足初始條件(8)的解T(t;T0)存在且唯一。 對于式(7),由引理2Ito公式有 根據Ito微分計算規則,結合式(7),可得 (9) 積分式(9),并由初始條件(8)得 即 (10) 由式(10)可知,當T0>0時,恒有T(t;T0)>0,?t>0。 證畢。 T(t;T0)的期望反映了全球平均溫度的變化情況,其方差往往反映了溫度的波動。下面分別討論T(t;T0)的期望ET(t;T0)和方差VarT(t;T0)。 定理2 對于方程(7),有 ET(t;T0)=T0exp(rt),VarT(t;T0)= 證明方程(7)的解可寫為如下的隨機過程 對該隨機過程兩邊取期望,并結合引理3,有 等式兩邊對t求導可得 因此,對方程進一步求解可得 ET(t;T0)=T0exp(rt) 對T2(t;T0)由引理2Ito公式可得 dT2(t;T0)=(2r+α2)T2(t;T0)dt+2αT2(t;T0)dB(t) 其隨機過程為 等式兩邊取期望有 同理,等式兩邊對t求導,再進一步求解可得 證畢。 考慮方程(7)的參數估計。令lnT(t;T0)=x(t),記lnT0=x0。 則式(9)變為 (11) 對于Δt>0,以等距時間點列0,Δt,2Δt,…,nΔt。 離散化方程(11)可得 其中εi(i=1,2,…,n)是獨立且均服從標準正態分布。因此由矩方程組(5)可設 解這個方程組,并且注意到有 因此 所以,有 其中,Ti表示第i(i=0,1,2,…,n)年的平均溫度。 由定理2知,當t→∞時,有ET(t;T0)→∞(t→∞), 即隨著時間的推移,溫度會無限增長。顯然這與實際情況不符,因為自然資源畢竟有限,有限的環境條件會對溫度的增長起到阻滯作用,全球溫度的上升會逐漸變慢并最終會在某個溫度附近波動。基于此,建立一個關于溫度的隨機阻滯增長模型 (12) 其中Tm表示環境容納量,即自然資源,環境條件所能容納的最高溫度。根據引理2Ito公式有 再將式(12)代入,可得 (13) 仍記方程(12)滿足初始條件(8)的解為T(t;T0),對式(13)積分有 (14) 證明當初始值T0>0時,根據式(14),顯然有T(t;T0)>0。再對式(12)進行變形,有 (15) 則根據式(15)有T(t;T0) 證畢。 受到自然資源等各種客觀因素的限制,環境容納量的值是一定的,因此,借助R語言軟件中的deSolve包對Tm值進行估計,得到的結果為Tm=16.6224。對于方程(12)其Euler顯示迭代公式為 T(tn)=T(tn-1)+rT(tn-1)· αT(tn-1)[B(tn)-B(tn-1)] 由此可以模擬得到含隨機擾動項和不含隨機擾動項的模型擬合效果(圖2)。 圖 2 模型擬合效果 圖2中,直線表示不含隨機擾動項的阻滯增長模型,曲線表示模型(11)所擬合的效果。由于內稟增長率r值較小,因此,不含隨機擾動項的模型增長緩慢,在短期內圖形呈直線增長。顯然本文擬合的隨機阻滯增長模型更加合理。 進一步,對全球平均溫度進行了一定量的預測(圖3)。 圖 3 全球平均溫度預測 由圖3可以看出,在未來的很長一段時間內,全球平均溫度還是會呈現緩慢、波動的上漲趨勢。通過數據擬合得到的環境容納量為16.6224℃,也就是說最終在上升了2℃左右溫度將趨于穩定。溫度首次達到這一數值的時間t=2395,也就是再過375年左右全球平均溫度才會趨于穩定,達到環境容納值。 全球平均氣溫上升2℃也是我們人類可以承受的極限。2015年11月30日在巴黎舉辦的聯合國氣候大會指出:如果到21世紀末,人類還無法將全球平均氣溫增幅控制在2℃以下,會導致海平面上升、森林火災、水資源缺乏等一系列災難性后果。這個結果也與2014年聯合國政府間氣候變化專門委員會(IPCC)在德國柏林發布的第五次評估報告的研究結果相同。這說明本文擬合出的環境容納量的值是準確有效的。
1.3 隨機微分方程的數值計算
1.4 矩法估計

2 模型的建立與分析
2.1 氣溫變化的隨機微分方程模型
2.2 氣溫變化的隨機指數增長模型









2.3 全球氣溫的隨機阻滯增長模型



3 數值模擬

