時(shí)間分?jǐn)?shù)階擴(kuò)散方程微分階數(shù)與擴(kuò)散系數(shù)聯(lián)合反演的貝葉斯方法

劉艷杰,王迎美,李功勝

(山東理工大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,山東 淄博 255049)

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

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

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

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

1 正問題與差分格式

1.1 正問題

考慮一個(gè)時(shí)間分?jǐn)?shù)階擴(kuò)散方程的初邊值問題：

(1)

(2)

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

1.2 差分格式

(3)

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

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

j=0時(shí)，

(4)

j>0時(shí)，

(5)

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

(6)

其中

(7)

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

其中：

當(dāng)i,l=1,2,…,M-1時(shí)，

(8)

當(dāng)i,l=2,…,M-2時(shí)，

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

(9)

關(guān)于差分格式(6)的穩(wěn)定性與收斂性，有下述結(jié)論：

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

2 反問題與反演模型

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

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

(10)

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

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

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

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

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

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

(11)

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

Y=F(θ)+ε。

(12)

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

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

3 貝葉斯推理與MCMC算法

3.1 貝葉斯推理

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

(13)

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

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

(14)

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

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

第2步：根據(jù)待估參數(shù)和測量數(shù)據(jù)的關(guān)系，提出一個(gè)似然函數(shù)p(Y|θ)；

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

第4步：在后驗(yàn)狀態(tài)空間進(jìn)行抽樣，建立后驗(yàn)分布的近似分布，取其統(tǒng)計(jì)量即為待估參數(shù)的反演值。

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

(15)

及

(16)

進(jìn)一步假設(shè)微分階數(shù)與擴(kuò)散系數(shù)相互獨(dú)立，則有：

(17)

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

(18)

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

(19)

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

3.2 MCMC方法

目前，MCMC方法已成為貝葉斯推理中后驗(yàn)抽樣的標(biāo)準(zhǔn)方法。在MCMC方法應(yīng)用中，有兩個(gè)較為重要的算法：Metropolis-Hasting算法(簡稱M-H算法)和Gibbs算法。本研究采用M-H算法對后驗(yàn)狀態(tài)空間進(jìn)行抽樣，算法步驟如下：

1) 在模型參數(shù)先驗(yàn)范圍內(nèi)，設(shè)定初始值θi(i=1)；

2) 利用當(dāng)前參數(shù)計(jì)算正問題，獲得該參數(shù)對應(yīng)的條件概率密度p(θi|Y)；

3) 根據(jù)參數(shù)的建議分布,在當(dāng)前參數(shù)狀態(tài)下提取新的測試參數(shù)θ*，并計(jì)算該參數(shù)對應(yīng)的條件概率密度p(θ*|Y)；

5) 重復(fù)步驟3)和4)，直至達(dá)到迭代次數(shù)。

4 數(shù)值反演

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

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

4.1 n=100,似然函數(shù)標(biāo)準(zhǔn)差的影響

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

圖1 反演誤差與似然函數(shù)的標(biāo)準(zhǔn)差Fig.1 Inversion error and standard deviation of likelihood function

(20)

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

分析圖1發(fā)現(xiàn)，標(biāo)準(zhǔn)差取8×10-5～1.5×10-3時(shí)，反演參數(shù)的均值誤差在較小的范圍內(nèi)波動，反演結(jié)果較好。取σ=1×10-4進(jìn)行反演計(jì)算，圖2(a)、圖2(b)分別繪制了反演值隨迭代次數(shù)的變化曲線與反演參數(shù)的后驗(yàn)直方圖。

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

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

表1 σ=1×10-4時(shí)的反演統(tǒng)計(jì)結(jié)果Tab.1 Inversion statistical results,when σ=1×10-4

4.2 σ=8×10-5,附加數(shù)據(jù)向量維數(shù)的影響

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

從表2的計(jì)算結(jié)果看出，附加數(shù)據(jù)向量的維數(shù)對反演算法的影響較小。隨著附加數(shù)據(jù)的減少，反演誤差與運(yùn)算時(shí)間的變化幅度均不大。同時(shí)表明，附加數(shù)據(jù)的個(gè)數(shù)不宜取得過多，這表現(xiàn)出與最佳攝動量算法等梯度型反演方法相似的特點(diǎn)。

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

5 結(jié)論

本研究應(yīng)用Bayesian-MCMC方法給出了時(shí)間分?jǐn)?shù)階擴(kuò)散方程微分階數(shù)與擴(kuò)散系數(shù)的統(tǒng)計(jì)反演，并探討了似然函數(shù)標(biāo)準(zhǔn)差與附加數(shù)據(jù)向量維數(shù)的不同選取對反演算法的影響。模擬計(jì)算結(jié)果表明，貝葉斯反演方法不依賴誤差泛函的梯度計(jì)算和未知量的初值選取，對于相對復(fù)雜模型的多參數(shù)反演是有效的。反演計(jì)算中似然函數(shù)的方差既不能太小也不能太大。當(dāng)方差落在一個(gè)可行區(qū)間時(shí)，反演結(jié)果較好。此外，附加數(shù)據(jù)向量的維數(shù)對反演算法幾乎沒有影響，但附加數(shù)據(jù)的個(gè)數(shù)不宜過多。