羅心悅,陳 瑜
(1. 復(fù)旦大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,上海 200433; 2. 上海財經(jīng)大學(xué) 數(shù)學(xué)學(xué)院,上海 200433)
植被斑圖(Vegetation pattern)是世界上諸多干旱半干旱地區(qū)常見的景觀特征,通常形成于年降雨量50 mm到750 mm的地區(qū),受降水、地質(zhì)地形、土壤和植被的空間異質(zhì)性及人類活動的影響,呈條帶狀、圓環(huán)狀、斑點狀等結(jié)構(gòu)。植被斑圖不僅能夠維持干旱半干旱地區(qū)能量與物質(zhì)的循環(huán)過程,還能防止風(fēng)蝕荒漠化。植被格局反映了干旱半干旱區(qū)生態(tài)系統(tǒng)的生態(tài)水文過程,為沙漠化等生態(tài)系統(tǒng)轉(zhuǎn)型提供預(yù)警指標(biāo)。因此,研究植被斑圖的組成結(jié)構(gòu)、形成機制以及演替過程,正確認(rèn)識生態(tài)系統(tǒng)的關(guān)鍵因子及各因子之間的相互關(guān)系,對于揭示干旱半干旱區(qū)生態(tài)系統(tǒng)變化的關(guān)鍵過程具有重要意義。這將有助于預(yù)測氣候變化對荒漠地區(qū)生態(tài)系統(tǒng)的影響,為地區(qū)生態(tài)環(huán)境建設(shè)形成綠色發(fā)展模式提供重要的理論基礎(chǔ)和技術(shù)支撐。
在理論分析方面,非線性科學(xué)的進步帶動了斑圖動力學(xué)的發(fā)展。反應(yīng)擴散方程是非線性系統(tǒng)研究中的一個重要分支,在化學(xué)實驗控制、生物識別、期權(quán)定價、影像分析等眾多前沿科學(xué)領(lǐng)域有重要應(yīng)用。反應(yīng)擴散模型能夠模擬反應(yīng)擴散過程中豐富的空間動態(tài)模式,因而其衍生出的各類問題一直為數(shù)學(xué)家和物理學(xué)家重點關(guān)注。學(xué)者們探究了生態(tài)水文過程、資源的空間再分配、植物個體間的空間相互作用,以及氣候和管理條件等因素對植被格局的影響[1]。即使在地形和土壤沒有表現(xiàn)出任何異質(zhì)性的情況下,仍可以觀察到植被格局的形成,這說明植被格局是內(nèi)在植被動力學(xué)的結(jié)果[2]。1999年,Klausmeier[3]提出刻畫斜坡上植被格局的數(shù)學(xué)模型,基于該模型對斜坡上呈條紋狀和點狀斑圖模式的成因進行解釋,并結(jié)合數(shù)值模擬得到條紋狀和點狀的空間斑圖。
目前對植被斑圖模型的數(shù)值研究主要集中在兩方面: 一是利用微分方程理論對模型的解進行定性研究,以及利用數(shù)值方法獲得模型近似解析解或高精度數(shù)值解;二是在已知模型并且獲得部分觀測數(shù)據(jù)后,確定最優(yōu)模型參數(shù),即參數(shù)反演。前者屬于正問題研究,后者在反問題研究范疇內(nèi)。觀測信息的不準(zhǔn)確和模型自身的耦合、非線性等特性,都為反問題的求解增加難度。因此,基于有限的觀測對模型參數(shù)進行準(zhǔn)確、穩(wěn)定地反演,識別出植被格局的演化特征,具有重大的現(xiàn)實意義。
考慮廣義Klausmeier-Gray-Scott模型。如圖1所示,模型刻畫帶有土壤水?dāng)U散反饋調(diào)節(jié)機制的平原地區(qū)植被斑圖的演化過程:
圖1 模型示意圖Fig.1 Model diagram
式中符號、參數(shù)的解釋說明如下:w、b為2維空間(x,y)隨時間t變化的水量和植被量;Δ是拉普拉斯算子,DwΔw、DbΔb分別表示水和植被的擴散,參數(shù)Dw、Db為各自的擴散系數(shù);A刻畫降水;Lw刻畫水汽蒸發(fā),參數(shù)L為蒸發(fā)水平;Rwb2表示支持植物生長的水份供給,也表示植物的生長,參數(shù)R為供給水平;Mb刻畫植被的消亡,參數(shù)M為植被的死亡率;Ω是邊界?Ω光滑的有界區(qū)域;初始條件φ(x,y)和ψ(x,y)是非負(fù)函數(shù);v表示水在負(fù)方向的下坡徑流量。
不失一般性[4],對參數(shù)空間進行簡化,令L=R=1,即有
(1)
本文采用了如下所示的1階半隱式后向歐拉差分格式(First-order semi-implicit backward Euler difference scheme)[5]:
在MATLAB程序中,空間步長取1,時間步長取0.05,水的擴散系數(shù)Dw取10,植被的擴散系數(shù)Db取1,總時長T取30。圖2為2維區(qū)域Ω=(0,400)×(0,400)內(nèi)終時植被的空間分布圖,A分別取0.48、0.52、0.45。我們可以看到當(dāng)降雨量A接近閾值(A=0.50)時,緩慢的植物擴散和快速的水分?jǐn)U散可以演化出一種臨界的植被狀態(tài)——植被集中于小區(qū)域。其中,藍(lán)色區(qū)域為裸土(b=0),紅色區(qū)域為高植被密度區(qū)域。當(dāng)降雨量增加時,不同初始狀態(tài)的植被分布演變出圖3所示的斑點狀、條帶狀、迷宮狀等具有自組織(Self-organized)特征的模式。
圖2 降水量處于閾值時聚集型植被分布Fig.2 Distribution of aggregated vegetation when precipitation is at threshold
圖3 不同形態(tài)的的植被模式Fig.3 Different vegetation patterns
圖3中,初始條件設(shè)定為
w(x,0)=0.6+0.15cos(2π(0.5x×y))+0.3sin(2π(0.5x+0.4y×R)),b(x,0)=0.5+0.12cos(2π(0.5x×y))+0.2sin(2π(0.4x+0.5y×R))。
式中:R是服從高斯分布N(0,1)的隨機數(shù)矩陣,且令初始矩陣中為負(fù)數(shù)的元素為零。再選取參數(shù)Dw=10,Db=1,A=0.9,M=0.45。令初始條件為
w(x,0)=0.2+R,b(x,0)=0.3+R。
圖4展示了該設(shè)定下植被斑圖的演變過程。從圖4中可以看到,在演變初期,系統(tǒng)變化較大,具有一定的振蕩,隨著時間推移,系統(tǒng)趨于平穩(wěn),斑圖形狀基本穩(wěn)定。
圖4 特定參數(shù)取值下植被格局的時空演變Fig.4 Spatiotemporal evolution of vegetation pattern under specific parameters
為探究不同參數(shù)對系統(tǒng)最終狀態(tài)的影響,本文對系統(tǒng)主要參數(shù)A,M進行了靈敏度分析,結(jié)果如圖5所示。縱坐標(biāo)為各格點的植被密度b的總和btotal,代表觀測區(qū)域的植被密度大小。
圖5 主要參數(shù)對整體密度的影響Fig.5 Influence of main parameters on overall density
從圖5可以看到,降水量、植物死亡率的變化能夠引起區(qū)域內(nèi)植被密度的大幅波動。因此,對A、M的精確識別對于植被系統(tǒng)的預(yù)測與控制具有重要的意義。
我們分別記降水量A和植物死亡率M兩個參數(shù)為
c1=A,c2=M。
根據(jù)觀測條件的不同,參數(shù)識別問題可以細(xì)分為終時觀測的參數(shù)反演問題與間斷觀測的參數(shù)反演問題[6]。顯然,反演參數(shù)A、M均具有上界,因此,我們規(guī)定參數(shù)的可行域:
Uad={(c1,c2)∈2:0 (2) 其中不等式右端的C1和C2是由有關(guān)圖靈系統(tǒng)的適定性分析所確定的。正問題陳述如下: 正問題對于可行域Uad中給定參數(shù)(c1,c2),計算滿足系統(tǒng)(1)的水濃度w(x,t)和植被密度b(x,t)。 (3) 其中:J被定義為 以反應(yīng)擴散系統(tǒng)(1)為約束條件。由γi加權(quán)的項度量了在T時的系統(tǒng)解與觀測值之間的差異。由λi加權(quán)的項有效地限制了參數(shù)c1,c2的大小,這既滿足式(2)的要求,又對觀測噪聲有一定容許度。通過使用λi加權(quán)項進行反演的過程被稱為Tikhonov正則化。使用Tikhonov正則化一定程度上可以避免偏微分方程相關(guān)反問題可能出現(xiàn)的不適定性[7]。對于正則化參數(shù)的選擇有諸多方法,例如L-曲線方法(L-curve method)。正則化參數(shù)的最佳選擇的討論可以參考文獻[8],本文不做贅述。通過選取正則化參數(shù),我們可以平衡方程的擬合與參數(shù)的大小控制?,F(xiàn)在,我們給出反問題的陳述: 根據(jù)最優(yōu)控制理論,本文引入了伴隨系統(tǒng)(Adjoint system): (4) 式中:p(x,t)(簡記為p),q(x,t)(簡記為q)為伴隨變量,且 由于伴隨方程在時間上是倒向的,因此需要最終條件而不是初始條件。最終條件為 此外,我們由伴隨系統(tǒng)得到關(guān)于伴隨變量的最優(yōu)控制的顯式表征,即最優(yōu)條件: (5) 狀態(tài)方程(State Equation, SE)(1)和伴隨方程(Adjoint Equation, AE)(4)以及最優(yōu)條件(5)稱為最優(yōu)系統(tǒng)(Optimality system)??梢宰C明映射(c1,c2)(w,b)是Gteaux可導(dǎo)的,并且其導(dǎo)數(shù)滿足 3.1.1 變步長梯度算法 為求解參數(shù)反演問題,本文首先考慮用文獻[10]中提到的變步長梯度算法(Variable step gradient algorithm)獲得最優(yōu)解和最優(yōu)參數(shù)的近似。算法對狀態(tài)方程、伴隨方程和終時T處的目標(biāo)函數(shù)進行了離散處理。具體地,算法首先對參數(shù)c1(0),c2(0)進行初始化設(shè)定,并且預(yù)設(shè)迭代步長λ以及用于檢驗成本函數(shù)收斂性的誤差容限ε。在每次迭代中,算法求解w(k),b(k)(k≥1)的非線性反應(yīng)擴散系統(tǒng),并儲存損失函數(shù)J(c1(k),c2(k))。與此同時,算法還計算了伴隨變量的值p(k),q(k),以及損失函數(shù)c1(k),c2(k)的導(dǎo)數(shù)dJ(c1(k),c2(k))/d(c1(k),c2(k))。如果成本函數(shù)減小,依照適當(dāng)?shù)牟介Lλ沿該方向繼續(xù)計算;如果成本函數(shù)沒減小,則減小步長λ。在確定好新步長后,通過標(biāo)準(zhǔn)的梯度算法對參數(shù)進行更新: 重復(fù)這一步驟,直至損失函數(shù)的相對誤差小于容限,即 式中:JN(k)指第k次迭代中時間節(jié)點N處的損失函數(shù)。但考慮到上述迭代算法可能出現(xiàn)停滯的狀況,我們設(shè)計異常處理: 如果步長小于誤差容限,即λ<ε,則退出算法。離散最優(yōu)控制算法的每次迭代都需要狀態(tài)方程的數(shù)值解到終時T,伴隨方程的數(shù)值解從終時T開始向前求解,梯度更新。該最優(yōu)化問題的標(biāo)準(zhǔn)算法計算成本巨大。但由于目標(biāo)函數(shù)的最優(yōu)解是系統(tǒng)的平穩(wěn)解(Stationary solutions),上述標(biāo)準(zhǔn)算法可以改進為: 選擇與目標(biāo)函數(shù)相等的狀態(tài)方程的初始值,取最終時間T僅等于2個時間步長τ并尋找對應(yīng)于平穩(wěn)初始數(shù)據(jù)的參數(shù)c1(k)和c2(k)。改進算法的其他步驟與標(biāo)準(zhǔn)算法步驟基本相同。改進算法的邏輯很容易理解,如果參數(shù)c1(k)和c2(k)是最優(yōu)的,則2個時間步2τ后的損失函數(shù)為零,算法停止;如果參數(shù)c1(k)和c2(k)是次優(yōu)的,則初始數(shù)據(jù)在2個時間步長2τ上演化,然后通過損失函數(shù)衡量系統(tǒng)解和目標(biāo)方程之間的差異,通過調(diào)整參數(shù)優(yōu)化變步長梯度算法。 3.1.2 數(shù)值結(jié)果 在數(shù)值實驗中,本文發(fā)現(xiàn)標(biāo)準(zhǔn)梯度算法的運算時間對初值的選取具有較強的依賴性,迭代過程容易陷入停滯,無法收斂。然而改進梯度算法收斂于幾乎所有初始值,并且能在較短時間內(nèi)精確估計出用于生成觀測值的參數(shù)值。表1顯示了標(biāo)準(zhǔn)梯度算法(參數(shù)用Asd,Msd表示)和改進梯度算法(參數(shù)用Amd,Mmd表示)的數(shù)值結(jié)果,后者具有更顯著的準(zhǔn)確率。 表1 參數(shù)反演結(jié)果Tab.1 Parameter inversion results 為驗證算法的收斂性,本文繪制了損失函數(shù)fLoss隨算法迭代的變化圖像,如圖6所示,k迭代表示算法的迭代次數(shù)。對于改進梯度算法,損失函數(shù)在最初的幾次迭代中快速下降,在此之后維持在較低水平;對于標(biāo)準(zhǔn)梯度算法,由于在每次梯度更新中算法都以恒定步長計算多步,所以損失函數(shù)整體趨勢在下降,但存在一定的振蕩。 圖6 算法的殘差-步數(shù)曲線對比Fig.6 Comparison of residual-iteration curve of two algorithms s取0對應(yīng)不含噪音的觀測值函數(shù)。圖7(見第330頁)展示了含1%水平噪音的觀測圖像以及根據(jù)改進梯度算法重構(gòu)出的終時植被斑圖。從圖中可以看到,重構(gòu)算法能夠較好地還原出原始觀測中植被濃度較高的區(qū)域,并保持植被分布的基本形狀。 圖7 終時植被斑圖的觀測及重構(gòu)Fig.7 Observation and reconstruction of terminal vegetation pattern 綜上,改進梯度算法能夠有效地從含噪終時觀測中得到參數(shù)估計。然而在現(xiàn)實中,植被演化的終時觀測是難以獲得的,更多的情形是,可觀測數(shù)據(jù)來源于遙感探測器,研究者們需要從離散的觀測數(shù)據(jù)中分析出系統(tǒng)演化特征,然后對未來發(fā)展趨勢進行預(yù)測。使用確定性算法,相當(dāng)于對每一間斷觀測進行最優(yōu)參數(shù)求解,沒有很好地利用觀測的歷史信息,甚至可能出現(xiàn)相鄰兩次求解結(jié)果懸殊較大的情況,因此需要更科學(xué)的方法為植被防護提供支持。 3.2.1 數(shù)據(jù)同化及貝葉斯理論基礎(chǔ) 數(shù)據(jù)同化(Data assimilation)是基于觀測信息結(jié)合數(shù)學(xué)物理模型模擬動態(tài)變化的一種方法,應(yīng)用于復(fù)雜系統(tǒng)的建模和對發(fā)展趨勢的預(yù)測,其目的是同化觀測數(shù)據(jù),為預(yù)測提供盡可能準(zhǔn)確的可信區(qū)間[11]。數(shù)據(jù)同化基于概率對模型狀態(tài)變量和參數(shù)進行最優(yōu)估計,為模型的規(guī)律探尋和預(yù)測提供了一條有效的途徑。在這個意義下,貝葉斯定理(Bayesian theory)是數(shù)據(jù)同化方法的理論基礎(chǔ)。 貝葉斯理論的主要思想是根據(jù)已有的觀測信息修正原有判斷,從而得到關(guān)于模型的后驗概率分布。對于參數(shù)反演,研究問題具有三大要素,分別是數(shù)據(jù)(觀測值)、關(guān)鍵參數(shù)和數(shù)學(xué)模型?;谄⒎址匠痰姆汉治隹蚣?數(shù)據(jù)和參數(shù)可以認(rèn)為是Banach空間中的元素,即假設(shè)U和Y是Banach空間,y∈Y是數(shù)據(jù),p∈U是參數(shù)。一般來說,p代表參數(shù)向量。當(dāng)參數(shù)值p確定時,可以對模型進行求解,計算出與該參數(shù)對應(yīng)的模型解。這個過程被稱為正問題的求解,而參數(shù)反演問題是一個反問題。正問題的求解可以表示為如下映射: G:U→Y。 由于G將參數(shù)映射到測量空間(Measurement space),因此G稱為觀測算子(Observation operator)。如果對于每一個參數(shù)組,模型存在唯一解,G在空間U中是良定義的;倘若G為多值函數(shù),則G在U中不是良定義的。 實際上,數(shù)據(jù)中常存在著一定水平的噪聲,可以通過概率分布來進行建模。設(shè)噪聲η服從概率分布0,進而觀測值可以表示為 y=G(p)+η。 (6) 式(6)表示了數(shù)據(jù)y的概率分布與參數(shù)p和噪聲η之間的關(guān)系,參數(shù)反演問題的三大要素串聯(lián)了起來。從建模的角度來看,只要數(shù)據(jù)y位于Banach空間中,式(6)具有良定義的性質(zhì)。在Banach空間的假設(shè)下,我們能夠得到參數(shù)反演問題的嚴(yán)格證明,并且結(jié)果適用于解的數(shù)值近似[12]。 在貝葉斯理論框架下,系統(tǒng)(1)的參數(shù)反演問題又可以陳述如下: 令μy為(p|y)的概率密度函數(shù)(Probability density function),(y|p)的概率密度函數(shù)由ρ(y-G(p))給出,其中ρ是噪聲分布0的概率密度函數(shù)。用μ0表示先驗分布0的概率密度函數(shù)。根據(jù)貝葉斯定理,我們得到 μy(p)∝ρ(y-G(p))μ0(p)。 貝葉斯定理表明,假設(shè)勢Φ:U×Y→是在乘積測度0×0可度量的,并假設(shè) 由于觀測值y是已知的,即Φ不依賴于y,所以我們可以用Φ(p)表示Φ(p;y)。用B(0,r)?U表示圓心在原點和半徑r的球,用 ‖ · ‖Σ表示帶有權(quán)重Σ的歐幾里得范數(shù)。利用正問題的性質(zhì),以下命題保證了參數(shù)反演問題的適定性: 后驗分布存在,且連續(xù)依賴于數(shù)據(jù)y。 3.2.2 數(shù)值結(jié)果 觀測數(shù)據(jù)的生成與上一節(jié)類似,初始條件均為空間均勻穩(wěn)態(tài)附加一定水平的擾動,參數(shù)取值為A=0.58、M=0.45、T=30。實驗設(shè)定10-5為穩(wěn)態(tài)閾值,即當(dāng)系統(tǒng)離散時間導(dǎo)數(shù)的L2范數(shù)小于閾值時,可以認(rèn)為系統(tǒng)達(dá)到了穩(wěn)態(tài)。實驗記錄了終時T以及T/10、T/5、T/3、2T/3時的數(shù)據(jù)。為了驗證系統(tǒng)解的穩(wěn)定性,實驗還繼續(xù)求解t=2T時刻的系統(tǒng)解,并檢查T時刻和2T時刻解之間的差值是否低于閾值。圖8表明,T時刻之后系統(tǒng)已達(dá)穩(wěn)態(tài)。 圖8 系統(tǒng)解的穩(wěn)定性Fig.8 Stability of system solution 實驗將標(biāo)準(zhǔn)偏差為5%的高斯噪聲添加到原始觀測數(shù)據(jù)(見圖9(a))中,所得結(jié)果如圖9(b)(第332頁)所示。 圖9 添加噪音后的觀測結(jié)果Fig.9 Observations after adding noise 參數(shù)的先驗信息是由[0.4,0.7]×[0.5,0.8]區(qū)域上的正態(tài)分布來建模的。在已知部分先驗信息之后,具體地,我們采用馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo, MCMC)抽樣算法,進行Metropolis Hastings采樣[13],計算模型參數(shù)的后驗分布。 針對正問題的求解耗時長的問題,本文在實現(xiàn)過程中引入了自適應(yīng)以及參數(shù)遷移的思想: 每次采樣模擬記錄下擬合表現(xiàn)最好的幾對參數(shù),并在下一次模擬中,讓模擬點盡可能與之相近,以期更快速地逼近最優(yōu)參數(shù)組。圖10展示了次數(shù)為100的自適應(yīng)MCMC采樣的樣本分布,圖11展示了加入隨機項后最后10次的采樣分布,k模擬為模擬次數(shù)。 圖10 自適應(yīng)MCMC采樣的樣本分布Fig.10 Sample distribution of adaptive MCMC sampling 圖11 加入隨機項后的采樣分布Fig.11 Sampling distribution after adding random items 在圖12中可以看到參數(shù)對(A,M)的后驗分布的95%置信域。仔細(xì)觀察置信域如何集中在精確值附近,我們會發(fā)現(xiàn)這不同于最初的認(rèn)識——在區(qū)域[0.4,0.7]×[0.5,0.8]上服從正態(tài)分布。數(shù)值實驗結(jié)果表明,95%置信域被包含在[0.587,0.616]×[0.438,0.458]范圍內(nèi),并且置信域在A方向的長度大約是M方向的1.5倍。此外,圖中信息還表明,可信域中心點距離真實的最優(yōu)值十分相近,這一定程度上驗證了算法的有效性。 圖12 參數(shù)A和M的后驗分布的95%置信域Fig.12 95% confidence region of posterior distribution for parameters A and M 圖13給出了基于反演所得參數(shù)值對系統(tǒng)進行重構(gòu)的結(jié)果。這些數(shù)據(jù)圖表明,本節(jié)算法可以很好地重構(gòu)植被斑圖的演化過程,預(yù)測結(jié)果與實際相符。 圖13 重構(gòu)植被斑圖的結(jié)果(上: 真實值,下: 重構(gòu)值)Fig.13 Results of reconstructing vegetation patch map 本文基于廣義Klausmiere-Gray-Scott方程,建立了帶有土壤水?dāng)U散反饋調(diào)節(jié)機制的平地植被斑圖的演化模型,分別對具有終時觀測和間斷觀測的參數(shù)識別問題,應(yīng)用了不同的反演算法,并在數(shù)值實驗中驗證了算法的有效性。植被呈現(xiàn)斑圖狀的地區(qū)的生態(tài)模式具有不穩(wěn)定性,植被斑圖的生成與演化受降水量、植被死亡率等自然因素的制約影響,所以關(guān)鍵參數(shù)識別對于植被防護、地理探測具有重要的現(xiàn)實意義。對于已知終時觀測的參數(shù)識別問題,變步長梯度算法能夠有效合理地估計出參數(shù),但算法具有一定的局限性: 不能用來估計與非平穩(wěn)目標(biāo)函數(shù)相關(guān)的參數(shù);如果損失函數(shù)中的目標(biāo)函數(shù)沒有接近反應(yīng)擴散系統(tǒng)的解,則算法可能無法收斂,產(chǎn)生不準(zhǔn)確的模型參數(shù)估計。而現(xiàn)實中,終時觀測難以獲得,觀測往往具有一定的離散特性,需要研究者從中分析系統(tǒng)演化特點,預(yù)判系統(tǒng)演化趨勢。貝葉斯理論為已知一定先驗信息的參數(shù)反演問題提供了途徑。本文使用馬爾可夫鏈蒙特卡羅(MCMC)抽樣方法來確定產(chǎn)生與數(shù)據(jù)相似模式的反應(yīng)擴散系統(tǒng)參數(shù)的分布,進行了若干次數(shù)值模擬,找到了參數(shù)的可信域。 由于實際區(qū)域的生態(tài)組成和觀測條件都較復(fù)雜,本文做了適當(dāng)?shù)暮喕?該模型有待進一步地被細(xì)化。植被斑圖的演化需要由耦合的非線性擴散反應(yīng)方程來刻畫,再之由于反問題自身的非線性,這為參數(shù)識別的提出與理論分析增加了一定的困難。如何給定適當(dāng)?shù)母郊訔l件保證所提出反問題的適定性,如何建立這類反問題的條件穩(wěn)定性,這都是值得進一步深入探討的問題。3.1 終時觀測的參數(shù)識別問題
3.2 間斷觀測的參數(shù)識別問題
4 結(jié) 語