時間分數階擴散方程微分階數與擴散系數聯合反演的貝葉斯方法

劉艷杰,王迎美,李功勝

(山東理工大學 數學與統計學院,山東 淄博 255049)

近二十年來，分數階微積分及反常擴散模型在物理學、力學、材料科學、環境科學與水文地質學等領域得到了廣泛的應用[1-4]。分數階反常擴散方程相關反問題研究受到越來越多的關注。對于時間分數階擴散方程，時間微分階數反映了反常擴散的遺傳特性和整體性，擴散系數主要反映了介質屬性和擴散速率，兩類參數在反常擴散模型模擬中起著非常重要的作用。實際的擴散問題中，微分階數與擴散系數往往都是未知的，需要依據擴散基本模型和附加信息加以識別確定，這就形成一類關于微分階數與擴散系數的雙參數聯合反演問題。對于時間分數階擴散方程單純確定微分階數或者確定擴散系數的反問題有一些研究[5-7]，但對于同時確定微分階數與擴散系數的反問題研究相對較少。文獻[8]證明了在脈沖初值條件下，應用單邊測量數據確定時間微分階數與空間依賴擴散系數反問題的唯一性；文獻[9]在一類逼近空間中，應用Levenberg-Marquart算法給出這類反問題的數值反演；文獻[10]考慮光滑初值條件下這類雙參數反問題的唯一性，并應用一種最佳攝動量正則化算法給出擾動數據下的數值反演。2015—2016年，文獻[11-12]對于一般時間分數階擴散方程同時確定多個微分階數與模型參數的反問題，應用解析方法證明了多參數聯合反演的唯一性。

上述研究或是給出了這類參數聯合反演問題的唯一性，或是應用梯度型算法及正則化策略進行了數值反演模擬，均是反常擴散相關反問題研究的主要方面。從模擬計算的角度看，傳統的反演方法大都基于誤差泛函的極小化，通過正則化策略、梯度計算與迭代方法獲得反演解的點估計，不僅算法依賴于梯度計算和初值選取，而且忽略了由于模型近似和數據擾動等因素造成的反演解的不確定性。近年來，隨著計算技術的進步和解決實際問題的需要，不確定性量化與統計計算方法研究受到關注，其中貝葉斯推斷方法在數學物理反問題研究中得到了廣泛應用。

貝葉斯方法主要通過融合先驗信息與樣本知識得到后驗分布，進而獲得研究對象的統計特征。由于后驗分布往往沒有顯性的表達式，借助馬爾可夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)抽樣技術的Bayesian-MCMC方法成為實現參數統計特征反演的主要途徑。Wang等[13-15]較早探討了參數反演問題的貝葉斯方法；Yan等[16]應用Bayesian-MCMC方法研究了一維熱傳導方程時間依賴Robin系數的統計特征反演；朱嵩等[17-18]研究了水文地質與水環境模型源項系數反演的貝葉斯算法；Yan等[19]通過構建替代模型提出貝葉斯反演的L1-SCM算法，并給出了二維時間分數階擴散方程的微分階數與點源位置參數聯合反演的算例；Jia等[20]應用貝葉斯方法研究了基于變指標Besov先驗的函數重建問題；Fan等[21]應用貝葉斯方法研究了一類反常分形擴散模型的多參數反演問題；Hu等[22]提出了無窮維貝葉斯推斷的一種自適應預處理Crank-Nicolson MCMC算法；Yan等[23]給出了基于替代模型的貝葉斯反演算法的收斂性。最近，Iglesias等[24]應用貝葉斯方法研究了一類帶有移動邊界的樹脂傳輸模型的參數反演；Zhang等[25]對于時間-空間分數階擴散方程微分階數和源項聯合反演的貝葉斯方法進行了研究；Hoang等[26]應用貝葉斯方法研究了一個橢圓方程具有局部周期的兩尺度系數的反演重建問題；Yan等[27]基于混合多項式展開技術，提出了參數反問題貝葉斯求解的多保真加速算法。

雖然貝葉斯方法應用于數理方程反問題的求解已有不少研究，但多側重于后驗分布的抽樣算法及其加速收斂問題或者先驗信息與正則化策略的融合問題，對于似然函數及其參數的選取以及附加數據對算法的影響等問題研究相對較少。本研究主要探究一個時間分數階擴散方程中微分階數與擴散系數聯合反演的貝葉斯方法。這類參數聯合反演問題在一定條件下具有唯一性，文獻[10]應用最佳攝動量正則化算法實現了不同參數取值條件下的數值反演。本研究應用Bayesian-MCMC算法給出這一反問題的統計特征反演，并討論似然函數的方差與附加數據量的選取對反演算法的影響。

1 正問題與差分格式

1.1 正問題

考慮一個時間分數階擴散方程的初邊值問題：

(1)

(2)

當微分階數α、擴散系數D以及初始分布g(x)均已知且滿足相容條件時，式(1)是一個適定的正問題，其數值求解已有不少研究。為了內容的完整性，下面給出式(1)的差分解法。

1.2 差分格式

(3)

相應地，式(1)中的初邊值條件離散為：

因此，可得隱式差分格式如下：

j=0時，

(4)

j>0時，

(5)

應用矩陣表示，式(4)～(5)即為：

(6)

其中

(7)

系數矩陣A=(ail)(M-1)×(M-1)，

其中：

當i,l=1,2,…,M-1時，

(8)

當i,l=2,…,M-2時，

aii=1+2p;a11=aM-1,M-1=1+p。

(9)

關于差分格式(6)的穩定性與收斂性，有下述結論：

引理1[28]對于正問題(1)，差分求解格式(6)是無條件穩定的。

2 反問題與反演模型

對于正演模型(1)，微分階數α∈(0,1)表示反常擴散的時間記憶特征，擴散系數D表示介質的屬性及擴散速率，這兩個參數在實際問題中都是難以直接測量獲取的。基于(1)，并根據額外的一些觀測信息，形成一個確定微分階數與擴散系數的反問題。文中考慮附加信息為區域右邊界處的觀測數據：

u(1,t)=h(t),0

(10)

對于由式(1)和式(10)形成的雙參數反演問題，根據文獻[8,10]中的研究結果可以證明反演解具有唯一性。

命題1設微分階數α∈(0,1)，擴散系數D>0，則雙參數反問題(1)和(10)具有唯一解。

證明：分兩種情況證明。設初值為特殊的Delta函數，即g(x)=δ(x)。注意到擴散系數D為一個正數，則根據文獻[8]的方法，即知α與D可由邊值數據u(1,t)(0

若初值取為一般的光滑函數，此時對于一般的擴散系數，反演的唯一性需要兩個邊界處的觀測數據或者擴散系數的額外信息。不過，注意到擴散系數恒為正數，其在區間[0,1]上滿足對稱性條件，則依據文獻[10]中的證明方法，可知該反問題的解唯一。證畢。

為便于討論和書寫，記θ=(α,D)及I=(0,1)×(0,θ)，這里θ<+∞為給定的正數。對于任意給定的θ∈I，求解正問題得到的解記為u(θ)(x,t)。結合附加條件(10)，反問題首先可化為時間域t∈(0,1]上一個方程的求解

h(t)=u(θ)(1,t)。

(11)

考慮到實際觀測數據和計算數據的離散性，上式左右兩端對應離散后分別記為Y與F(θ)，則得到一個有限維的隨機反演模型

Y=F(θ)+ε。

(12)

其中：F(θ)是給定參數時的正問題計算值；θ=(α,D)是需要評估反演的模型參數；Y=h(t)表示觀測數據向量；ε表示反演模型的隨機誤差，包括數據測量誤差、正問題計算誤差及計算機舍入誤差等，一般服從于均值為0、方差為σ2的高斯分布。

以下主要應用貝葉斯方法求解模型(12)，獲得參數θ=(α,D)的分布值及其統計特征。

3 貝葉斯推理與MCMC算法

3.1 貝葉斯推理

對于反演模型(12)，記p(θ)是待估參數的先驗分布概率密度函數，p(Y|θ)是似然函數，表示給定參數取值時觀測數據分布的條件概率，則有貝葉斯公式

(13)

其中p(θ|Y)表示參數的后驗分布概率密度函數。這里的反演即是要求出在給定觀測數據Y的條件下使得后驗分布概率最大的參數分布。由于p(Y)是積分常量，具體計算中(13)式常寫為：

p(θ|Y)∝p(θ)p(Y|θ)。

(14)

貝葉斯反演是基于(14)式的一種參數估計方法，其實施步驟為：

第1步：依據待估參數的先驗信息，得到先驗概率密度函數p(θ)；

第2步：根據待估參數和測量數據的關系，提出一個似然函數p(Y|θ)；

第3步：根據(14)式計算得到待估參數的后驗概率密度函數p(θ|Y)；

第4步：在后驗狀態空間進行抽樣，建立后驗分布的近似分布，取其統計量即為待估參數的反演值。

根據先驗信息，獲得待估參數的先驗概率分布是實現貝葉斯算法的前提。歷史經驗或專家經驗是重要的參考，常用的先驗分布有高斯分布、均勻分布和馬爾可夫隨機場等。注意到這里的待估參數是落于區域I的正常數，假設均服從均勻分布，即α～U(0,1),D～U(0,5)，則有：

(15)

及

(16)

進一步假設微分階數與擴散系數相互獨立，則有：

(17)

似然函數表征了模型參數擬合觀測數據的程度，其值越大則擬合效果越好。在反演模型(12)中，假設誤差ε服從均值為0、方差為σ2的正態分布，則有：

(18)

其中n表示觀測數據向量Y的維數，即觀測數據的多少。綜上，根據貝葉斯公式即得：

(19)

通過(19)式即可計算出參數的后驗概率密度函數。然而，由于正演關系F(θ)較為復雜難以解析表達，使得后驗概率密度p(θ|Y)很難直觀表示出來。為了獲得參數估計值，需要采用特定的隨機抽樣方法去獲得待估參數的近似分布，再利用這個近似分布的統計量作為待估參數的估計值。常用的隨機抽樣方法是馬爾可夫鏈蒙特卡洛方法。

3.2 MCMC方法

目前，MCMC方法已成為貝葉斯推理中后驗抽樣的標準方法。在MCMC方法應用中，有兩個較為重要的算法：Metropolis-Hasting算法(簡稱M-H算法)和Gibbs算法。本研究采用M-H算法對后驗狀態空間進行抽樣，算法步驟如下：

1) 在模型參數先驗范圍內，設定初始值θi(i=1)；

2) 利用當前參數計算正問題，獲得該參數對應的條件概率密度p(θi|Y)；

3) 根據參數的建議分布,在當前參數狀態下提取新的測試參數θ*，并計算該參數對應的條件概率密度p(θ*|Y)；

5) 重復步驟3)和4)，直至達到迭代次數。

4 數值反演

本節應用上述Bayesian-MCMC算法對參數反問題進行數值反演。不失一般性，取微分階數真值為α=0.8，擴散系數真值D=1.5，即θ=(0.8,1.5)。利用該參數真值及差分格式(6)求解正問題，得到附加數據向量Y，另取待估參數的建議分布為均勻分布：q(θ*|θi)=U(θi-Δθ,θi+Δθ)，且設步長為先驗范圍的5%。

對于本研究的雙參數反演問題，有兩個參數的選取值得關注，一是向量Y的維數,即附加數據的個數n；二是似然函數的標準差σ。此外，若無特殊說明，反演中的參數初始值均在先驗限定范圍內隨機產生，反演迭代總次數取為10 000次。首先考察使用較多附加數據情況下，似然函數的標準差對反演算法的影響。

4.1 n=100,似然函數標準差的影響

設附加數據向量Y的維數dim(Y)=100，即n=100。又設似然函數的標準差σ∈[2×10-5,5×10-3]，對于不同的標準差分別進行反演計算得到微分階數與擴散系數的后驗均值。圖1繪制了聯合反演誤差error隨似然函數的標準差σ的變化曲線，這里聯合反演誤差表示為：

圖1 反演誤差與似然函數的標準差Fig.1 Inversion error and standard deviation of likelihood function

(20)

其中E(·)表示后驗均值。

分析圖1發現，標準差取8×10-5～1.5×10-3時，反演參數的均值誤差在較小的范圍內波動，反演結果較好。取σ=1×10-4進行反演計算，圖2(a)、圖2(b)分別繪制了反演值隨迭代次數的變化曲線與反演參數的后驗直方圖。

圖2 σ=1×10-4時反演迭代曲線與后驗直方圖Fig.2 Inversion iterative curve and posterior histogram,when σ=1×10-4

從圖2(a)可以看出，馬爾可夫鏈在經歷了2 000次迭代后即達到收斂狀態，微分階數與擴散系數均收斂到真值左右。從圖2(b)看出，參數反演值分別出現在0.8和1.5左右的概率最大，這與參數真值吻合。此外，采用后8 000次的反演數據統計計算兩個參數的中位數和均值，統計結果列于表1。從表1看出，中位數和均值反演值均收斂于參數真值，微分階數反演的最大誤差不超過0.204 3%，擴散系數反演的最大誤差不超過0.778 8%。

表1 σ=1×10-4時的反演統計結果Tab.1 Inversion statistical results,when σ=1×10-4

4.2 σ=8×10-5,附加數據向量維數的影響

取定σ=8×10-5，其他參數及符號表示同上，考察附加數據向量維數變化對反演算法的影響。為了均衡計算過程的隨機性，計算結果取5次反演的平均值。表2給出了附加數據向量維數n=100,n=40,n=20與n=10時的平均反演結果，其中Err表示平均反演誤差，Tcpu表示平均的CPU運行時間。

從表2的計算結果看出，附加數據向量的維數對反演算法的影響較小。隨著附加數據的減少，反演誤差與運算時間的變化幅度均不大。同時表明，附加數據的個數不宜取得過多，這表現出與最佳攝動量算法等梯度型反演方法相似的特點。

表2 σ=8×10-5,附加數據向量維數與反演結果Tab.2 Inversion results with dimension of additional data vector,when σ=8×10-5

5 結論

本研究應用Bayesian-MCMC方法給出了時間分數階擴散方程微分階數與擴散系數的統計反演，并探討了似然函數標準差與附加數據向量維數的不同選取對反演算法的影響。模擬計算結果表明，貝葉斯反演方法不依賴誤差泛函的梯度計算和未知量的初值選取，對于相對復雜模型的多參數反演是有效的。反演計算中似然函數的方差既不能太小也不能太大。當方差落在一個可行區間時，反演結果較好。此外，附加數據向量的維數對反演算法幾乎沒有影響，但附加數據的個數不宜過多。