Lorenz系統參數估計方法研究

章翠蓮，李維德*，朱高峰

(1.蘭州大學數學與統計學院, 甘肅 蘭州 730000; 2.蘭州大學資源環境學院, 甘肅 蘭州 730000)

Lorenz系統參數估計方法研究

章翠蓮1，李維德1*，朱高峰2

(1.蘭州大學數學與統計學院, 甘肅 蘭州 730000; 2.蘭州大學資源環境學院, 甘肅 蘭州 730000)

馬爾科夫鏈蒙特卡洛方法(簡稱MCMC)可用于復雜系統的不確定性估計和參數估計.基于貝葉斯理論，運用幾種改進的MCMC方法：自適應Metropolis 算法(AM)、延遲拒絕自適應Metropolis算法(DRAM)、差分進化馬爾科夫鏈算法(DE-MC)對具有復雜的動態性質的Lorenz 混沌系統未知參數進行了探討性的估計.根據未知參數的后驗概率密度似然函數，利用MATLAB 仿真，選取樣本點出現頻率高的區間作為目標分布區域，并通過縮小先驗分布范圍來計算參數的估計值.分析比較這三個算法模擬的結果，得出如下結論：合適的目標分布區域在參數估計中很關鍵；在搜索樣本點上，DRAM 算法的遍歷性最高；DE-MC 算法可使Lorenz 系統獲得較為精確的參數向量，更適用于復雜系統未知參數的估計.

AM；DRAM；差分進化馬爾科夫鏈；Lorenz 混沌系統；參數估計

1963年，洛倫茲提出了第一個表現奇異吸引子的動力學系統，即Lorenz系統[1].該系統表現出非線性動力學系統的復雜形式——具有混沌動態性質.近年來，混沌系統在保密通信，人工智能以及生命科學等領域得到應用和發展.但由于混沌系統的復雜性，參數的不可觀測性，或者保密性，系統的參數是未知的.因此系統的參數估計在混沌控制以及同步領域上的作用是不可估量的[2].

通過觀測數據求模型的未知參數，是建模分析中的重要步驟，但由于觀測過程中存在系統誤差以及觀測信息的偶然誤差，不一定可求得合適的參數解，或者得出的參數解不唯一，因此本文采用基于貝葉斯理論下改進的MCMC方法對Lorenz系統未知參數估計進行探討性的研究.

MCMC的Metropolis-Hasting(M-H)算法，其思想是構造一個以目標分布π(x)為不變分布的馬爾科夫鏈,通過一個概率密度函數q(x,y) (稱為建議函數)來實現這一分布[3,4]. 該方法是否具有收斂性取決于是否選取合適的建議分布函數.假設當前狀態為x，則從建議函數q(x,.)產生一個候選點y，那么接受概率:

(1)

若α>u(u是[0,1]均勻分布的隨機數)，接受候選點y作為新的采樣點，否則復制原來的x作為新的采樣點.

當建議分布為對稱分布時，即q(x,y)=q(y,x).此時接受概率為：

(2)

近年來，許多研究人員在MCMC算法的運行效率和效果做出許多努力從而改進了方法.HeikkiHaario等提出了自適應Metropolis算法，簡稱AM.這個算法利用了MarkovChain的歷史信息自動調節提議函數的協方差矩陣，與隨機游走算法相比，提高了算法運行效率[5,6].文獻[2] 應用了該算法對Lorenz系統的參數進行了估計.在AM算法提出后，HeikkiHaario等又提出了自適應Metropolis采樣和延遲拒絕相結合的算法：延遲拒絕自適應Metropolis算法，簡稱DRAM.這個算法提高了馬爾科夫鏈的遍歷性和采樣點的接受率[7].TerBraak則通過遺傳算法思想，提出了一種差分進化馬爾科夫鏈算法，簡稱DE-MC.該算法多條鏈同時運行，其解決了MCMC的一個重要問題，即為跳躍分布提供了合適的范圍和方向[8].基于貝葉斯理論，本文主要將這三種改進的MCMC算法運用到對Lorenz系統的未知參數估計，對這三種方法進行探討性的討論，得出解決Lorenz系統未知參數估計的較優方法.

1 理論與算法

1.1 貝葉斯理論

貝葉斯方法提供了一個涵蓋模型不確定性的框架，在這個框架下，具有概率分布的貝葉斯方法的關鍵的特點在于描述了參數和模型的不確定性[9].

基于貝葉斯理論，我們認為未知參數向量(對Lorenz系統而言，其有3個未知參數，因此d=3),在對X進行初步推斷時，根據相關未知參數的信息 ，我們認為其先驗分布為p(θ).在對系統進行觀測采樣后，得到觀測信息(包含N個樣本點).因此，由貝葉斯估計，系統參數的后驗分布，觀測信息和先驗分布之間的關系[2]：

(3)

考慮到Z是已知的信息，上式可以寫成：

p(θ|Z)∝p(Z|θ)P(θ)

(4)

我們知道，貝葉斯推斷的核心在于使用似然函數去分析參數的不確定性.因此對p(Z|), 在實際中一般用似然函數表示，其值越大表明擬合程度越好[10].

1.2AM算法

AM是基于Metropolis算法中的隨機游走算法(RWM)[11]改進而來的,該算法關鍵在于，先構造具有高斯分布的建議函數，利用馬爾科夫鏈先前全部樣本信息計算建議分布的協方差矩陣，自適應地向目標函數逼近[12]. 利用該點，通過后驗分布獲得的信息使得建議分布得到更新.在第i步中，Haario等提出了以當前點作為均值和具有協方差矩陣Ci的多元正態分布.該協方差Ci的定義為：

(5)

其中，sd是縮放因子，a1>0 是一個很小的數，它們分別取為2.382/d與10-6.在理論上要保證Ci非奇異，Id是d維單位矩陣.C0是初始協方差矩陣.Cov(θ0,...,θi-1)表示θ0,...,θi-1的協方差矩陣. 若初始協方差過大也會使AM算法自適應性變差；反之，初始協方差過小降低AM算法的遍歷性. 為非自適應長度，其值越大自適應性越慢[12].

AM算法流程如下：

a)設定i=1，對所求的參數變量初始化；

b)利用(5)式計算協方差矩陣Ci；

c)產生候選點θ*～N(θi-1,Ci)；

d)根據(1)式，計算接受概率α=min{1,p(θ*|Z)/p(θi-1|Z)}；

e)產生滿足均勻分布的隨機數u～U[0,1]；

f)若α>u,則θi=θ*，否則θi=θi-1；

g)重復b)～f)，直到完成產生指定的迭代數量的樣本.

1.3DRAM算法

DRAM是基于延遲拒絕算法[13]和AM算法相結合而得來的算法[7].

DR算法是具有不同提議函數和轉移核的MCMC方法.Peskun和Tierney已證明出有限維的狀態空間和普遍的狀態空間里DR算法通過降低停留在當前候選點的概率提高MCMC方法的運行效率[14,15].其主要思想在于：在進行M-H抽樣時，若當前的候選點被拒絕，不是保持當前狀態不變，而是從當前位置利用改變的提議函數再次移動，通過對第2次移動的候選點的接受概率計算來保持馬爾科夫鏈的可逆性，其本質上是通過對提議函數局部自適應調整來提高算法的效率[12].

在算法進程中，假定當前位置為x，用π(·)表示該點的目標分布函數，根據多元正態分布函數q0(x,·)產生新的候選點y，接受率(類似式(1))為：

(6)

若在點x被拒絕，則根據當前位置和被拒的候選點再次移動，從新的建議函數q1(x,y,·)產生新的候選點z,相應的接受概率為[16]：

(7)

若是第二級的候選點仍被拒絕，則可在被拒的候選點z上再進行第三級移動.

DRAM與AM算法的流程大致相同，不同于AM算法的是在第f)步.將第f)步改動:若α>u, 則θi=θ*，再返回到第b)步繼續迭代；若α

1.4DE-MC算法

差分進化馬爾科夫鏈是在具有MCMC模擬的實參數空間上利用Metropolis準則與差分進化算法[17]結合而得來的MCMC算法[8].該算法克服了MCMC方法中建議分布需要選取合適的范圍和方向以及相應的計算效率問題.該算法有兩個參數：縮放因子γ和平行鏈的數目N可根據實際情況定義[18].

DE-MC算法的基本過程[8]：

θ*=θi+γ(θa-θb)+ε

(8)

其中，θ*是建議選取的參數集，θi是當前的參數集，θa,θb是從θi不含外的參數集隨機選取出來的，γ為縮放因子.ε～N(0,b)d,在這里b很小.建議的參數集接受或者拒絕取決于Metropolis比率,具體可參考(1) 式[16].

上面所列舉的三種MCMC方法，AM和DRAM算法都具有自適應性，DE-MC是遺傳算法和傳統的Metropolis方法相結合的算法.每一種算法各有其相關的優勢，因此下面將這三種算法運用到Lorenz模型進行未知參數的估計，進行探討性的討論.

2 Lorenz混沌系統參數估計分析

Lorenz系統可由如下的常微分方程組表示：

(9)

我們知道當系統(9)中a,b,r分別為10,28,8/3時，系統表現為混沌現象.Lorenz描述的非線性耗散系統具有初始條件敏感性.在用數值解法求解非線性微分問題時，數值積分方法和編程現實中的誤差累積導致模擬的軌道與原來的軌道偏離.即不同的初始條件，不同的算法和編程細節，使得計算軌道是不可重復的.因此，系統的長期行為是不可預測的.本文選取的時間長度為10s.首先對該系統進行數據采樣：先讓Lorenz系統進行演化,以混沌系統中某一點(1.500，2.230，-1.270)作為初始點，并記為：(x1(0),x2(0),x3(0)).設置時間步長h為0.01s,采用四階定步長Runge-Kutta算法來求解常微分方程，得出離散的數值解，我們得到Lorenz系統在10s內的1000組離散數據： (x1(ti),x2(ti),x3(ti))(i=0,1,...,1000).

它們可作為觀測數據集Z, 由于這些數據在操作過程中必有觀測誤差，因此可多模擬幾次，取每個步長狀態的平均值作為觀測數據集來減少誤差的干擾.

現在已知觀測數據集，分別采用AM,DRAM和DE-MC估計出Lorenz系統的未知參數集θ=(a,b,r).由貝葉斯公式求出參數的聯合后驗概率密度函數.令模型的初始狀態量為：(x1(0),x2(0),x3(0))，每隔10個步長共取10個樣本點，即抽取0h，10h，...，90h時刻的狀態量(x1(ti),x2(ti),x3(ti))(i=1,2,...,10)作為這次的觀測數據集Z.

根據貝葉斯公式(3),未知參數的后驗分布：

p(a,b,r|Z)=k·p(Z|a,b,r)p(a,b,r)

(10)

其中，k是一個大于0的常數.

在貝葉斯推斷中，首先設定未知參數的先驗分布p(a,b,r).假設這三個參數均滿足獨立均勻分布，用U(a),U(b),U(r)分別表示三個參數的均勻分布，則

p(a,b,r)=U(a)U(b)U(r)

(11)

相比較其他分布，均勻分布是最簡單的分布，其限定了參數變化的界限，先驗分布的界限很小時，目標區域則相對變小，這樣有利于提高參數的精確度.通常情況下熟知歷史經驗更能得到精確的先驗分布.

對于式(10)中的p(Z|a,b,r),前面提到觀測集是經過多次測量取平均值而得到的.根據普適似然不確定性分析方法GLUE的基本原理,用似然函數來比較模擬值和實測值的擬合程度，本文選用基于殘差和的似然函數[19,20]：

(12)

本文取n=10,d=3,M=1.對向量進行更新后，通過用四階定步長Runge-Kutta算法解Lorenz系統微分方程，時間和步長不變，每隔10h取一個數據，則得到10個數據.第i個數據記為:(y1(ti),y2(ti),y3(ti))(i=1,2,...,10).

由(10)-(12)可得，該系統未知參數的后驗概率密度似然函數為：

(13)

在模擬之前，本文先假設三個參數的目標分布區間均為[0,30],即a,b,r的聯合先驗分布均為：

(14)

后驗概率密度似然函數確定后，則可以用AM,DRAM,DE-MC算法對這三個未知參數進行隨機抽樣.對于AM和DRAM算法,非自適應段長度=200.DRAM算法中，λ=0.1，m=3.

3 仿真結果及分析

考慮到三個參數的目標分布都較大，迭代次數先取為2000，參數初始值則取相對應目標分布中的隨機值. 在每個算法進行抽樣時，每次都對未知參數向量進行更新，通過計算該隨機參數向量的似然后驗密度函數值與上一步后驗密度函數值的比值，決定是否將該向量作為下一步的值并進行更新.

在各參數的目標分布區間皆為[0,30]時，各個算法模擬出來各參數的馬爾科夫鏈還未收斂(圖1)，仿真的結果不理想.根據圖1中基于殘差平方和似然函數值較大的部分對應的參數向量，即a=9.9529,b=28.2056,r=2.7923時，式(12)中的較大，說明我們估計的參數向量在這個向量的附近.但由于MCMC算法具有隨機性，三種算法模擬出來三個參數的直方圖都不理想(圖2)， 縮小各個參數的目標分布區域再進行Matlab模擬可期望得出更精確的參數向量.重新設置參數的目標分布：a～U[9.9529-2,9.9529+2],b～U[28.2056-2,30],r～U[2.7923-2，2.7923+2]，其他條件不變.

圖1 目標分布區間為[0,30]下各算法仿真結果

注：T表示迭代次數，Q表示目標分布區域，L表示基于殘差和的似然函數值；(1)-(3)分別是用AM模擬得的參數a,b,r的馬爾科夫鏈；(5)-(7)分別是用$DRAM$模擬得的參數a,b,r的馬爾科夫鏈；(9)-(11)分別是用DE-MC模擬得參數a,b,r的馬爾科夫鏈；(4),(8),(12)分別是用AM,DRAM,DE-MC模擬得到的基于殘差和的似然函數值演化過程.

圖2 目標分布區間為[0,30]下各算法仿真的參數值方圖

注：P表示頻率；(1)-(3)分別是用AM模擬出來參數a,b,r的直方圖；(4)-(6)分別是用DRAM模擬出來參數a,b,r的直方圖；(7)-(9)分別是用DE-MC模擬出來參數a,b,r的直方圖.

圖3 第一次調整目標分布后各算法仿真結果

在新的參數目標分布區域下，由三個算法模擬出各參數的馬爾科夫鏈以及計算出式(12)的值的演化過程(圖3)，我們可以看出DE-MC算法模擬出來的參數r的馬爾科夫鏈已收斂，其他兩個參數在小范圍里波動，其基于殘差的平方和的似然函數值大體比其他兩個算法的大.在迭代次數1700-1800范圍內，式(12)的值最大，此時參數向量為(9.9875,28.1047,2.6508)，AM和DRAM算法模擬各參數的馬爾科夫鏈皆無收斂的趨勢.從遍歷性來看，DRAM的遍歷性最好(候選點接受率約為0.5)，AM算法次之(候選點接受率約為0.35)，DE-MC的遍歷性最差.結合三個算法模擬出各個參數的直方圖(圖4)，無法從AM算法模擬出樣本點的直方圖減小目標分布區域；DRAM算法模擬出來的各參數樣本點直方圖類似正態分布，參數r的最為明顯；從DE-MC算法模擬出的各參數樣本點直方圖可以看出樣本點皆分布集中，因此可根據其縮小目標分布區域.考慮DE-MC算法中的式(12)似然函數較大值所對應的參數向量，可將參數向量的目標分布設置為：(9.9875±0.5,28.1047±0.5,2.6508 ±0.2).

圖4 第一次目標分布調整后各算法仿真的參數值方圖

利用調整后的目標分布區域得出三個算法模擬出來的各參數的馬爾科夫鏈以及由式(12)得到的似然函數值演化圖(圖5).從圖5可看出，DRAM算法與其他兩個算法相比，其遍歷性最高，樣本點的接受概率約為0.54,AM算法樣本點的接受概率越為0.14，但是這兩個算法最后未收斂.從后驗概率密度似然函數來看，DE-MC算法收斂，并且其似然函數值最大.選取似然函數值最大的部分，即用DE-MC算法模擬出的馬爾科夫鏈收斂的那部分，其對應的參數向量約為：(9.9949, 28.0042, 2.6673), 即約為9.9949，約為28.0042，約為2.6673. 將其帶入(9)式，經過數值模擬，由(12)式可得出 約為1006，即其與觀測集的殘差平方和約為0.001. 由此可見，DE-MC算法更適用于Lorenz混沌系統的未知參數估計.

圖5 第二次目標分布調整后各算法仿真結果

圖6 狀態變量估計值與真實值在15-25s內的軌跡對比圖

Matlab中解一階微分方程組初值問題常見的方法為ode45函數，該函數實現了變步長四階五級的Runge-Kutta-Felhberg算法，采用變步長算法解微分方程[21].將原來的參數向量和通過DE-MC算法求解出來的參數向量代入Lorenz模型的狀態方程(9)，利用MATLAB中ode45函數數值解法求解.當時間長度為10s，初始狀態相同，發現各狀態變量真實值和估計值在10s內的振幅以及頻率相似，相位大體相同，軌跡幾乎一模一樣.其他條件不變，延長時間長度，各狀態變量真實值和估計值在20s后開始不在一條軌道上，估計值漸漸脫離原來的軌道(圖6).因此，Lorenz系統的長期行為是不可預測的.在選取觀測樣本數據時，系統不宜長時間演化，時間一般不超過10s. 在適當的目標分布區域內，DE-MC算法適用于短時間演化的復雜系統的參數估計，其最后模擬出的參數向量是有意義的.

4 結論與討論

本文基于貝葉斯理論的框架，利用MCMC方法中的AM,DRAM和DE-MC算法對Lorenz模型未知參數向量進行采樣，經過兩次目標分布區域的調整，分析比較各參數的馬爾科夫鏈，得到了系統的參數估計值.本文研究的主要結論總結如下：

(1)當未知參數的目標分布區域較大時，這三個算法模擬出各參數的馬爾科夫鏈皆不收斂，基于殘差和的似然函數值不高，效果不是很理想.參照各算法模擬出的似然函數值演化過程圖以及各參數樣本點的直方圖，選取樣本點出現頻率高的區間作為新的目標分布區域，及時調整先驗分布.參數的先驗分布越準確，即目標分布區域越小，模擬出來的參數向量越精確.

(2)本文應用了MCMC的三種算法對Lorenz混沌模型未知參數進行了估計.與文獻[2]仿真的結果不同的是，從運行過程來看，在相同的迭代次數下，AM算法運行的時間最少；由于DE-MC算法是多條鏈同時運行，模擬過程最緩慢.從樣本點的搜索區域來看，DRAM算法搜索范圍大，遍歷性最好；而DE-MC算法的遍歷性較差.從模擬出來的結果來看，經過兩次目標分布區域的調整，在三種算法中，DE-MC算法模擬出的馬爾科夫鏈收斂，并且基于殘差和的似然函數值較大，說明模擬出來的參數向量與真實的參數向量相差甚小.

(3)從參數估計的效果看，在適當的目標分布區域內，對于短時間演化的Lorenz混沌系統，DE-MC算法能較快地模擬出該參數向量.對一般的非線性以及動態的模型，DE-MC算法同樣能估計模型的參數.

綜上所述，AM,DRAM和DE-MC算法皆是通過構造馬爾科夫鏈進行隨機模擬，在對Lorenz混沌系統進行未知參數估計的過程中，實際就是在先驗分布界定的目標區域對參數向量進行隨機且較好抽取樣本點的過程，也就是構造馬爾科夫鏈的過程.通過Matlab仿真模擬，DE-MC算法能較好地對Lorenz混沌系統未知參數進行估計，因此該算法能更好地應用在復雜系統的參數估計上.

[1]劉福才, 梁曉明,Hénon. 混沌系統的廣義預測控制與同步快速算法[J].物理學報,2005, (10):4584-4589.

[2]曹小群,宋君強.基于MCMC方法的Lorenz混沌系統的參數估計[J].國防科技大學學報,2010,(2):68-72.

[3]ChristropheA,NandoDF.AnintroductiontoMCMCformachinelearning[J].MachineLearning,2003,50,5-43.

[4]陳平，徐若曦.Metropolis-Hastings自適應算法及其應用[J].系統工程理論與實踐,2008，(1):100-108.

[5]HaarioH,SaksmanE,TamminenJ.AnadaptiveMetropolisalgoritnm[J].Bernoulli,2001,(2):223-242.

[6]SaksmanE,ViholaM.OntheergodicityoftheadaptiveMetropolisalgorithmonunboundeddomain[J].TheAnnalsofAppliedProbability,2010,(6):2178-2203.

[7]HaarioH,LaineM,MiraA,etal.DRAM:EfficientadaptiveMCMC[J].Statistics&Computing, 2006,(4):339-354.

[8]BraakCJ.AMarkovChainMonteCarloversionofthegeneticalgorithmdifferentialevolution:EasyBayesiancomputingforrealparameterspaces[J].Statistics&Computing, 2006, (16):239-249.

[9]LombardiMJ.Bayesianinferenceforstableditributions:ArandomwalkMCMCapproach[J].ComputationalStatistics&DataAnalysis,2007,(5):2688-2700.

[10]MarshallL,NottD,SharmaA.AcomparativestudyofMarkovchainMonteCarlomethodsforconceptualrainfall-runoffmodeling[J].WaterResourcesResearch, 2004, (40):183-188.

[11]KuczeraG,ParentE.MonteCarloassessmentofparameteruncertaintyinconceptualcatchmentmodels:TheMetropolisalgorithm[J].JournalofHydrology, 1998, (1-4):69-85.

[12]郝燕玲,單志明.基于DRAM算法的α穩定分布參數估計[J].華中科技大學學報，2011,(10):73-78.

[13]TierneyL,MiraA.SomeadaptiveMonteCarlomethodsforBayesianinference[J].StatisticsinMedicine, 1998, (17-18):2507-2515.

[14]PeskunPH.OptimumMonte-CarlosamplingusingMarkovchains[J].BiometrikaTrust, 1973,(3):607-612.

[15]TierneyL.AnoteonMetropolis-Hastingskernelsforgeneralstatespaces[J].AnnalsofAppliedProbability,1998, (1):1-9.

[16]SmithTJ,MarshallLA.Bayesianmethodsinhydrologicmodeling:AstudyofrecentadvancementsinMarkovchainMonteCarlotechniques[J].WaterResourcesResearch, 2008, (12):67-76.

[17]StornR,PriceK.Differentialevolution-Asimpleandefficientheuristicforglobaloptimizationovercontinuousspaces[J].JournalofGlobalOptimization, 1997,(4):341-359.

[18]BraakCJFT,VrugtJA.DifferentialevolutionMarkovchainwithsnookerupdaterandfewerchains[J].Statistics&Computing, 2008,(4):435-446.

[19]BevenKJBinleyAM.Thefutureofdistributedmodels:Modelscalibrationanduncertaintyprediction[J].HydrologicProcess,1992,(2):279-298.

[20]FreerJK,BevenKJ.Bayesianestimationofuncertaintyinrunoffpredictionandthevalueofdata:AnapplicationoftheGLUEapproach[J].WaterResourcesResearch,1996,(7):2161-2173.

[21]薛定宇,陳陽泉.高等應用數學問題的MATLAB求解(第3版)[M].北京：清華大學出版社，2013.

(責任編校：晴川)

Study on Methods of Parameters Estimation of Lorenz System

ZHANG Cuilian1,LI Weide1*, ZHU Gaofeng2

(1.School of Mathematics and Statistics, Lanzhou University, Lanzhou Gansu 730000, China; 2. School of Resources and Environment, Lanzhou University, Lanzhou Gansu 730000, China)

Markov Chain Monte Carlo method (in short as MCMC) is useful to the uncertainty and parameter estimation for complicated system. Based on the Bayesian theory, the paper proposes an explorative study by utilizing several MCMC methods, such as Adaptive Metropolis (AM), Delay Rejection Adaptive Metropolis (DRAM) and Differential Evolution Markov Chain (DE-MC) algorithm, to estimate the parameters of Lorenz system which has complicated dynamical property. According to the posterior probability density likelihood function of unknown parameter, utilizing MATLAB simulation, selecting the region where sampling points appearing with high frequency as target distribution area, and narrowing the scope of the prior distribution, we calculate the estimation value of parameters. By analyzing and comparing the simulated results of the three algorithms mentioned above, we achieved these conclusions: suitable target distribution area is critical to estimate parameter; DRAM algorithm has high ergodicity in searching sampling points; DE-MC algorithm is more appropriate for complicated system with unknown parameter estimation, which makes Lorenz system obtain precise parameter vector.

AM; DRAM; differential evolution Markov Chain; Lorenz chaotic system; parameter estimation

2016-10-14

國家自然科學基金(批準號：41571016)資助項目.

章翠蓮(1991— )，女，廣西欽州人，蘭州大學數學與統計學院碩士生.研究方向：應用數學.

O212.1

A

1008-4681(2017)02-0005-06

*通訊作者：李維德(1967— )，男，甘肅蘭州人，蘭州大學教學與統計學院教授，博士.研究方向：數字生態學.