環(huán)境污染與自然資源耦合系統(tǒng)的動(dòng)力學(xué)模型分析

顧恩國，魯嘉珺

(中南民族大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)學(xué)院，武漢 430074)

環(huán)境污染與自然資源耦合系統(tǒng)的動(dòng)力學(xué)模型分析

顧恩國，魯嘉珺

(中南民族大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)學(xué)院，武漢 430074)

為研究環(huán)境自凈率和污染物排放速度對自然資源存量的影響，建立了自然資源與環(huán)境污染耦合系統(tǒng)的離散動(dòng)力學(xué)模型,討論了正平衡點(diǎn)的存在性、穩(wěn)定性及局部分叉，應(yīng)用中心流形理論探討了不動(dòng)點(diǎn)在發(fā)生flip分叉時(shí)的臨界穩(wěn)定性，并用數(shù)值模擬的方法驗(yàn)證了正平衡點(diǎn)的局部分叉.最后，用全局分析方法描述了環(huán)境自凈率對自然資源可持續(xù)利用的影響.

自然資源；污染物；正平衡點(diǎn)的穩(wěn)定性；flip 分叉；中心流形；可行吸引域

AbstractA discrete dynamic model of the coupled system of natural resources and environmental pollution was established, and the influence of environmental self-purification rate and pollutant emission rate on the stock of natural resources was analyzed. The existence, stability and local bifurcations of positive fixed points were investigated, the critical stability of fixed point in flip bifurcation was analyzed by using the center manifold theory, and numerical simulation was used to verify the local bifurcations of fixed point. At last, the global method was used to analyze the influence of environmental self-purification rate on the sustainable utilization of the natural resources.

Keywordsnatural resources; pollutants; stability of the positive equilibrium point; flip bifurcation; center manifold; feasible domain

隨著人口的增長和社會(huì)經(jīng)濟(jì)的發(fā)展，對生態(tài)系統(tǒng)造成的沖擊越來越強(qiáng).自然資源和服務(wù)的可持續(xù)性問題, 其重要性和緊迫性已經(jīng)引起了生態(tài)學(xué)家和經(jīng)濟(jì)學(xué)家的廣泛關(guān)注[1,2]，自然資源的可持續(xù)利用已經(jīng)成為非常重要的研究課題.在國內(nèi)外,關(guān)于自然資源可持續(xù)發(fā)展的研究很多,學(xué)者們建立了大量數(shù)學(xué)模型[3-6]，將自然資源與環(huán)境污染耦合起來建立動(dòng)力學(xué)模型并且用離散動(dòng)力學(xué)方法進(jìn)行研究的成果卻很少,對其演化進(jìn)行全局分析的更微乎其微.陳六君在文獻(xiàn)[7]中建立了自然資源環(huán)境系統(tǒng)的連續(xù)動(dòng)力學(xué)模型,并對模型所反映的突變性質(zhì)及其經(jīng)濟(jì)意義進(jìn)行了分析.本文在假設(shè)人類提取自然資源是自然資源儲(chǔ)量的線性函數(shù)的情況下，建立了自然資源儲(chǔ)量與污染物濃度相互作用的二維離散動(dòng)力學(xué)模型，用來研究環(huán)境自凈率和污染物排放速度對自然資源儲(chǔ)量的影響，進(jìn)而討論對自然資源可持續(xù)利用的影響.

1 模型的建立

本文假設(shè)自然資源增量在不考慮生物死亡和人類索取情況下滿足Logistic律：

N(t+1)-N(t)=rN(t)[1-N(t)/K],

其中N(t)為在時(shí)刻t自然資源存量，r為自然資源再生率；K為最大自然資源存量.假設(shè)污染物濃度增量與人類向環(huán)境系統(tǒng)輸入污染物的速度和環(huán)境自凈率有關(guān):D(t+1)-D(t)=p-cD(t)N(t),D(t)為在時(shí)刻t污染物濃度；p為人類向環(huán)境系統(tǒng)輸入污染物的速度;c為單位自然資源存量的自凈率.假設(shè)人類對自然資源的提取速度H(t)與自然資源存量呈線性關(guān)系，令H(t)=αN(t)，其中α為資源提取的速度.同時(shí)污染濃度會(huì)影響自然資源的再生能力，利用文[7]中的Hill函數(shù)，假設(shè)生物死亡函數(shù)：F(t)=fN(t)D(t)/(D(t)+D0)，f為致死率.根據(jù)上述假設(shè)和說明，建立系統(tǒng)模型如下：

(1)

2 正平衡點(diǎn)的存在性和穩(wěn)定性

二維動(dòng)力系統(tǒng)(1)可以寫成一個(gè)映射動(dòng)力系統(tǒng)：

(2)

根據(jù)平衡點(diǎn)的定義，系統(tǒng)(2)的平衡點(diǎn)就是下面非線性方程組的解：

(3)

由方程組(3)得：

(4)

(5)

聯(lián)立(4)和(5)式得到方程cK(r-μ-f)D2+[cK(r-μ)D0-μf]D-rpD0=0.

記a=cK(r-μ-f)≥0,b=cK(r-μ)D0-μf,Δ=b2+4arpD0>0.

關(guān)于正平衡點(diǎn)的存在性，可有下面的命題1.

由命題1知道系統(tǒng)(2)存在一個(gè)唯一的正不動(dòng)點(diǎn)，但這不足以得出自然資源可持續(xù)利用.為了得到自然資源可持續(xù)利用的條件，討論系統(tǒng)(2)的唯一平衡點(diǎn)的穩(wěn)定性.

系統(tǒng)(2)的Jacobian 矩陣是：

注意(4)和(5)式，Jacobian矩陣化簡為：

如果所有特征值在復(fù)平面上位于單位圓內(nèi)，則不動(dòng)點(diǎn)漸近穩(wěn)定，由Jury條件，若不動(dòng)點(diǎn)滿足條件：

則系統(tǒng)(2)正平衡點(diǎn)是漸近穩(wěn)定的.這3個(gè)條件對應(yīng)著特征值在復(fù)平面上以3種方式越過單位圓：若p1<0,至少有一個(gè)特征值大于1；若p2<0至少有一個(gè)特征值小于-1；若p3<0,有一對共軛復(fù)特征值位于單位圓外.它們等于0分別對應(yīng)著fold、flip、Neimark分叉.

綜上所述，可得到關(guān)于正不動(dòng)點(diǎn)穩(wěn)定的命題2.

取參數(shù)K=4,r=2,f=0.2,D0=0.5,μ=0.3,p=0.2,c=0.4,顯然r=2>0.5=f+μ, 正不動(dòng)點(diǎn)E1(3.3071,0.1512)漸近穩(wěn)定.圖1給出了參數(shù)為K=4,r=2,f=0.2,D0=0.5,μ=0.3,p=0.2時(shí)p1(c),p2(c),p3(c),00.5525，正不動(dòng)點(diǎn)E1通過flip分叉失去穩(wěn)定產(chǎn)生穩(wěn)定的2周期軌道，這與圖2關(guān)于c的一維分叉圖一致. 圖3(a)給出了關(guān)于c,p的二維分叉圖，黑色區(qū)域表示吸引子為不動(dòng)點(diǎn)對應(yīng)的參數(shù)范圍，白色區(qū)域?yàn)?周期吸引子對應(yīng)的參數(shù)范圍，花色區(qū)域表示吸引子周期大于20或混沌吸引子對應(yīng)的參數(shù)范圍，淺灰區(qū)域表示無窮吸引子對應(yīng)的參數(shù)范圍(表示資源爆炸). 我們發(fā)現(xiàn)并不是自凈率越高越好，當(dāng)自凈率超過某一值時(shí)，資源會(huì)爆炸. 圖3(b)給出了參數(shù)為K=4,r=2,f=0.2,D0=0.5,μ=0.3,c=0.57關(guān)于p的一維分叉圖,當(dāng)0.62

圖1 p1(c),p2(c),p3(c),0

圖2 關(guān)于c的一維分叉圖Fig.2 One dimensional bifurcation diagram with c

圖3 一維或二維分叉圖 Fig.3 One or two dimensional bifurcation diagram

3 產(chǎn)生flip分叉時(shí)正不動(dòng)點(diǎn)的穩(wěn)定性分析

當(dāng)參數(shù)集滿足p2=0時(shí)，在正不動(dòng)點(diǎn)E1處出現(xiàn)flip分叉.接下來通過中心流形理論確定正不動(dòng)點(diǎn)的穩(wěn)定性.取K=4,r=2,f=0.2,D0=0.5,μ=0.3,p=0.2,c=0.5525,此時(shí)參數(shù)集滿足p2=0, 正不動(dòng)點(diǎn)為E1(3.3285,0.1087).

首先，將系統(tǒng)(2)的正不動(dòng)點(diǎn)E1通過平移換轉(zhuǎn)換到原點(diǎn)O(0,0)，則原系統(tǒng)(2)化為：

所以系統(tǒng)(6)在E1(3.3285,0.1087)處的Taylor展式(忽略高階余項(xiàng))為：

(7)

將上述方程組寫成矩陣形式為：

(8)

其中系數(shù)矩陣為:

(9)

其中：

根據(jù)中心流形的定義，令：

v=h(u)=αu2+βu3+O(u4),

(10)

則h(Mu+f1(u,h(u)))-Nh(u)-g1(u,h(u))=0.

(11)

其中，M=-1，N=0.2238，f1(u,h(u))=0.4989u2+0.1408αu3+O(u4)，h(Mu+f1(u,h(u)))=h(-u+f1(u,h(u)))=αu2+(0.9978-β)u3+O(u4),

Nh(u)=0.2238(αu2+βu3+O(u4)),g1(u,h(u))=0.0338u2-0.4152αu3+O(u4),

將上式代入方程(11)得：

求得中心流形為：

h(u)=0.4355u2+0.5952u3+O(u4).

(12)

系統(tǒng)(9)在中心流形Mc的限制為:

f:u|→-u+f(u,h(u))=-u+0.4989u2-

0.0641u3+O(u4).

(13)

映射f在原點(diǎn)處的Schwarzian導(dǎo)數(shù)為:

(14)

因此h(u)為穩(wěn)定流形，從而正不動(dòng)點(diǎn)E1穩(wěn)定.

4 系統(tǒng)的全局分析

通過對正平衡態(tài)可行吸引域的性質(zhì)的分析，可以得出自然資源的自凈率在怎樣的條件下是維持資源可再生能力還是最終導(dǎo)致資源枯竭.本文我們主要采用數(shù)值模擬的方法來分析系統(tǒng)(2)的可行吸引域結(jié)構(gòu)， 根據(jù)系統(tǒng)(2)所代表的實(shí)際意義， 我們僅研究軌線完全在Ω={(N,D)|N>0,D∈R}中的N和D滿足的初始條件.若資源的存量是非正值，對自然環(huán)境資源的可持續(xù)利用將無法實(shí)現(xiàn)，通過數(shù)值試驗(yàn)，可以研究隨著自然資源自凈率變化的可行吸引域的變化情況.

圖4給出了自然資源可行吸引域的吸引子，選取的參數(shù)除p,c外與局部分析的值一樣，即K=4,r=2,f=0.2,D0=0.5,μ=0.3. 在圖4(a)-(c)中，黑色表示有界吸引子，白色區(qū)域表示有界吸引子的吸引域，深灰區(qū)域表示無窮吸引子的吸引域. 圖4(d)表明p=0.69,c=0.55時(shí)，不動(dòng)點(diǎn)與周期2軌道兩個(gè)有界吸引子共存，白色區(qū)域表示不動(dòng)點(diǎn)(空心點(diǎn))的吸引域，淺灰區(qū)域表示周期2軌道(實(shí)心點(diǎn))的吸引域.

圖4 可行吸引域和吸引子Fig.4 Feasible domains of attraction and attractors

5 結(jié)語

本文建立了自然資源與環(huán)境污染耦合的離散動(dòng)力學(xué)模型，對模型進(jìn)行了非線性動(dòng)力學(xué)分析，討論了正平衡點(diǎn)的存在性、穩(wěn)定性及局部分叉，并對正平衡態(tài)進(jìn)行了全局分析.研究發(fā)現(xiàn)：系統(tǒng)(2)具有混沌和共存吸引子復(fù)雜的動(dòng)力學(xué)行為，并且自然資源自凈率和人類向環(huán)境系統(tǒng)輸入污染物的速度對可再生資源的可持續(xù)利用具有顯著影響，當(dāng)它們都不太大時(shí)，系統(tǒng)(2)始終穩(wěn)定在正平衡點(diǎn)(見圖3(a))，并具有較大的可行吸引域的面積(見圖4(a))，因此針對不同的自然環(huán)境因素，要保證自然資源的可持續(xù)利用，需要結(jié)合自然資源的自凈率，控制人類向環(huán)境系統(tǒng)輸入污染物的速度，這對于自然資源管理決策者制定資源利用、保護(hù)政策具有重要的參考意義.

[1] Food and Agriculture Organization. The state of world fisheries and aquaculture[R]. Sofia: FAO, 2004.

[2] Garcia,Grainger J R.Gloom and doom？ The future of marine capture fisheries[J].Philosophical Transactions of the Royal Society,2005, B360:21-24.

[3] Clark C W. Mathematical bioeconomics: the optimal management of renewable resources[M].New York:Wilery,1976.

[4] Bischi G I,Kopel M,Szidarovszky F. Expectation-stock dynamics in multi-agent fisheries[J].Annals of Operations Research,2005,137:299-300.

[5] Kar T K,Matsuda H. A bioeconomic model of a single-species fishery with a marine reserve[J].Journal of Environmental Management,2008,86:171-180.

[6] Bischi G I, Szidarovszky F. Harvesting dynamics in protected and unprotected areas[J].Journal of Economic and Organization,2007,62:370-384.

[7] 陳六君，方?？? 自然資源環(huán)境系統(tǒng)的突變機(jī)制[J].復(fù)雜系統(tǒng)與復(fù)雜性科學(xué),2004,1(2):1-8.

DynamicModelAnalysisofCoupledSystemofEnvironmentalPollutionandNaturalResources

GuEnguo,LuJiajun

(College of Mathematics and Statistics, South-Central University for Nationalities, Wuhan 430074, China)

O192；F062.1

A

1672-4321(2017)03-0142-05

2017-06-20

顧恩國(1964-)，男,教授,研究方向：非線性動(dòng)力系統(tǒng)應(yīng)用，E-mail: guenguo@163.com

國家自然科學(xué)基金資助項(xiàng)目(61374085)；中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)基金資助項(xiàng)目(CZP17016)