時間分數階二維對流擴散方程多點源強的數值反演

李慧玲， 李功勝， 賈現正， 池光勝

(1.山東理工大學 理學院， 山東 淄博 255091；2.山東凱文科技職業學院 本科教育學院, 山東 濟南 250200)

分布式參數系統模型的建立對精確描述和有效控制許多物理和工程問題起關鍵作用.源項的估計或識別問題在環境科學和工業應用等領域顯得尤為重要.譬如，在環境水力學領域，如何尋找河流、城市水環境和湖泊的污染源；在化學反應過程、微波加熱過程和工業設計與制造領域，如何探測未知的發熱源等.

對拋物型方程的點源識別反問題，近年來國內外不少學者從不同的角度進行了研究.文獻[1]考慮了兩端無污染的流域中單個污染源的識別問題，在對污染源及測量點的先驗假設下，證明了源強識別的唯一性，并給出了識別算法.文獻[2-3]考慮了對流擴散方程的點源反演問題，利用某個時刻空間點處的測量數據將穩恒點源反演轉化為優化問題進行求解.文獻[4]將拋物型方程轉化為雙曲型方程后，證明了由內部測量數據反演多點源的唯一性.文獻[5]在第一類Dirichlet邊界條件下，利用熱變換方法將熱傳導方程多點源反問題轉化為等價的雙曲型方程反問題，然后通過分析該雙曲型方程反問題得到原反問題的唯一性和條件穩定性.文獻[6]研究了出流端為第二類邊界和零初始條件下單個污染源的識別問題，證明了源項反演的唯一性，并給出了一種局部穩定性和識別算法，文獻[7]討論了隨時間變化的單個污染源的反演問題及其應用.最近，文獻[8]研究了一維擴散方程中多個點源強度的識別問題，利用最佳攝動量正則化算法進行了有效的數值反演模擬.

本文將探討時間分數階二維對流擴散方程中確定多個點源的反演問題.對于這類分數階擴散方程點源識別反問題的研究尚未見有文獻報道.文中將首先給出正問題及其數值求解的差分格式，然后引入最佳攝動量正則化算法，并對多個點源強度進行精確數據和擾動數據條件下的數值反演.反演結果表明對于時間分數階二維對流擴散問題，應用最佳攝動量正則化算法可以實現對多個點源強度的數值確定.

1 正問題及其數值求解

考慮矩形域Ω=(0，l1)×(0,l2)上的時間分數階二維擴散方程

(1)

其中(x,y)∈Ω,t>0;D>0是擴散系數，v>0是沿著X軸方向的對流速度，g(x,y,t)是線性源項，而α∈(0,1)是時間分數階導數的階數.當α=1時,方程(1)即為通常的整數階二維對流擴散方程.對于帶非零源項的方程(1),給定零初邊值條件:

(2)

這樣,由方程(1)及初邊值條件(2),就構成一個二維時間分數階對流擴散的正問題. 下面給出這個正問題求解的一個有限差分格式.

1.1 正問題求解的差分格式

O(τ).

(3)

對于方程中的整數階導數項，按照通常的差分離散方法(參見文獻[11])，即有

(4)

(5)

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

當n=0時

(6)

當n>0時

(7)

1.2 數值試驗

表1 終值時刻(T=1)的精確解與數值解的比較

表2 給定空間點(x=y=0.5)處的精確解與數值解的比較

表3 不同微分階數對正問題求解的影響(x=y=0.5,T=1)

從表1、表2的計算結果可以看出,無論是在終值時刻,還是在給定的空間點處,正問題的精確解與數值解都吻合得較好.但從表3看出,時間微分階數對正問題求解具有一定的影響.微分階數越小,解誤差越小,當其接近于1時,解誤差逐步變大.

2 多點源項識別反問題與反演算法

2.1 反問題的提出

很多實際問題中,方程(1)中的源項是點源分布,且一般具有形式

(8)

其中：q為點源的個數；(xs,ys)為點源的位置；Qs為相應的點源強度(單位時間排放量)；δ為狄拉克函數.

假設已經知道點源的個數及其位置坐標，那么為了確定各個點源的強度值，需要補充關于污染物濃度分布的附加條件，并聯合正問題(1)-(2)形成一個源強識別反問題.本文給定t=T時的觀測值為附加數據,記u(xi,yj,T)=uT(xi,yj),對于每個固定的i,令j從1取到K,可定義附加數據向量VT：

VT

(9)

這樣，由附加數據(9)聯合正問題(1)-(2)構成了一個確定源強度Qs(s=1,2,…，q)的多點源強識別反問題.下面給出最佳攝動量正則化算法，并對上述源強識別反問題進行數值反演模擬.

2.2 最佳攝動量正則化算法

記Q=[Q1,Q2,…，Qq],利用上一節的差分方法求解正問題可得其解,并在t=T時刻賦值,記之為u(x,y,T;Q),稱為對應于輸入數據Q=[Q1,Q2，…，Qq]的計算輸出.

(10)

根據最佳攝動量算法(參見文獻[8-10, 12])，上述極小問題(10)的求解又轉化為對于給定的Qn,通過求解最佳攝動量δQn進而確定Qn+1的一種迭代算法：

Qn+1=Qn+δQn,n=0,1,2,…

(11)

且δQn是下述目標函數的極小值

F(δQn)=

(12)

將u(xi,yj,T;Qn+δQn)在Qn處作泰勒展開得到

u(xi,yj,T;Qn+δQn)=u(xi,yj,T;Qn)+

▽Qnu(xi,yj,T;Qn)·δQn+o(δQn),

略去高階項，則目標函數F(δQn)近似可得

u(xi,yj,T;Qn)-uT(xi,yj)]2+

μ‖δQn‖2

(13)

令ai,j=▽Qnu(xi,yj,T;Qn)，bi,j=uT(xi,yj)-u(xi,yj,T;Qn)其中

ai,j(s)=

λ為數值微分步長. 于是

μ‖δQn‖2

(14)

再根據最小二乘法的思想，求解minF(δQn)相當于求解規范方程

[ATA+μI]δQn=ATB

(15)

以下給出反演Q的算法步驟:

(1)給定未知量Q的初始猜測向量Q0和數值微分步長λ，求向量A,B.

(2)選取正則參數μ，求解方程(15)得到擾動量δQn=[ATA+μI]-1ATB.

(3)對于給定精度eps，判定是否滿足‖δQn‖≤eps.若是，則Qn即為所求，算法終止；否則，由(11)式得到Qn+1，再轉到步驟1繼續進行.

3 數值反演

本節應用最佳攝動量算法對源項識別反問題(1)-(2)及(9)進行數值反演。以下若無特殊說明，正問題計算中均取l1=l2=1，擴散系數D=0.001，平均流速v=0.001，終值時刻T=10，且M=10,N=100.

3.1 兩個點源的情形

設附加數據在(0.1, 0.1), (0.1, 0.2), (0.1, 0.3), (0.1, 0.4), (0.1, 0.5), (0.1, 0.6), (0.1, 0.7), (0.1, 0.8)以及(0.1, 0.9)等9個點處取得.假設在(0.2, 0.2), (0.8, 0.8)處有兩個點源，其源強分別為2和10，即源強真值為Q=[2,10].

(1)正則參數對反演結果的影響

取數值微分步長λ=0.1，初始迭代值Q0=[0,0]，反演計算結果列于表4，其中μ是正則參數，Qinv表示反演解，Err=‖Qinv-Q‖2表示反演解與源強真解的誤差，j表示迭代次數.

表4 正則參數對反演結果的影響

由表4可以看出正則參數越小反演結果越好，當μ=0時反演結果最好.當不用顯式正則化時(μ=0)，反演結果最好.以下若無說明，反演計算中都將取正則參數為零，且不再列出迭代次數.

(2)微分階數對反演結果的影響

取數值微分步長λ=0.1，初始迭代值Q0=[0,0]，計算結果見表5.

表5 微分階數對反演結果的影響

通過表5可以看出微分階數對反演結果的影響不大.

(3)初始迭代對反演結果的影響

取微分階數α=0.6，數值微分步長λ=0.1，計算結果見表6.

表6 初始迭代對反演結果的影響

通過表6可以看出初始迭代值對反演結果的影響不大.

(4)數值微分步長對反演結果的影響

取微分階數α=0.6，初始迭代值Q0=[0,0]，計算結果見表7.

表7 數值微分步長對反演結果的影響

通過表7可以看出數值微分步長對反演結果也幾乎沒有影響.

3.2 三個點源的情形

設在點(0.2,0.4),(0.4,0.6),(0.8,0.4)處有三個污染源，其源強分別為2，-5和10，即源強真值為Q=[2,-5,7].取D=1,v=0.1,Q0=[0,0,0].

(1)擾動水平對反演結果的影響

仍取微分階數α=0.6，數值微分步長λ=0.1，分別取擾動水平δ=10%,5%,1%進行反演計算,10次反演的平均值列于表8.

表8 擾動水平對反演結果的影響

通過表8可以看出隨著數據擾動水平的減小，反演結果精度越來越高.

(2)分數微分階數對反演結果的影響

設擾動水平δ=1%，考察分數微分階數對反演算法的影響，計算結果見表9.

表9 給定擾動水平下分數微分階數對反演結果的影響

由表9可以看出,對于給定的擾動數據,分數微分階數對反演結果的影響不大.

3.3 四個點源的情形

設在(0.2, 0.4), (0.4, 0.6), (0.6, 0.7)及(0.8, 0.8)處有4個點源，其源強分別為2，5，-7和10，即源強真值為Q=[2,5,-7,10].D=1,v=0.1,Q0=[0,0,0,0].

(1)附加數據的選取對反演結果的影響

取微分階數α=0.7，數值微分步長λ=0.1，附加數據選取x=0.1直線段上各點(0.1,yk)(k=1,…，K)處的觀測值.應用精確數據進行反演，對應于不同觀測點的計算結果見表10，其中y=[·]表示在x=0.1線上所取觀測點的縱坐標.

表10 附加數據的選取對反演結果的影響

通過表10可以看出，附加數據對反演結果具有一定的影響.若獲得附加數據的觀測點少于4個,則反演計算難以得到理想結果.但當選取4個以上觀測數據點時,反演結果相差并不大.

(2) 附加數據的擾動對反演結果的影響

仍在x=0.1直線段上選取觀測點,不妨設有9個觀測點,即x=0.1,y=[0.2,0.3,…，0.9], 考察附加數據的隨機擾動對反演算法的影響.十次反演計算的平均結果列于表11.

表11 附加數據擾動對反演結果的影響

通過表11可以看出,在四個點源的情形,數據擾動對反演算法影響較大,這說明此時的反演問題病態性較重.不過,隨著擾動水平的減小,反演誤差也變得越來越小.

4 結束語

本文探討了二維時間分數階對流擴散方程中多點源強度的識別反問題.基于正問題求解的差分格式,應用最佳攝動量算法進行了有效的數值反演.反演結果表明這類源項反問題的病態性并不強,特別當點源個數較少時,完全可以采用不帶顯式正則化項的最佳攝動量算法進行反演確定.對于四個點源的情形,只要給定四個附加數據點就可以實施模擬計算,但此時附加數據的擾動對反演算法影響較大,這也說明在點源個數多于四個的情形,反演問題的病態性增大,需要尋找更為有效的反演方法才能獲得更穩定的結果.

[1]王澤文,徐定華.流域點污染源識別的唯一性[J].寧夏大學學報：自然科學版,2006,27(2):124-129.

[2]閔濤,周孝德,張世梅,等.對流-擴散方程源項識別反問題的遺傳算法[J].水動力學研究與進展,2004,19(4):520-524.

[3]韓龍喜.河道一維污染源控制反問題[J].水科學進展,2001,12(1):39-44.

[4]支元洪,江茂澤,穆春來.拋物型多點熱源反演的唯一性[J].四川大學學報,2007,44(3):454-460.

[5]張小明,王澤文,王燕等.一類熱傳導方程多點源反演的唯一性和穩定性[J].華東理工大學學報(自然科學版) ,2009,32(2):189-193.

[6] El Badia A, Ha-Duong T, Hamdi A. Identification of a point source in a linear advection-dispersion-reaction equations: application to a pollution source problem[J].Inverse Problems, 2005, 21: 1 121-1 136.

[7] Hamdi A. Identification of a time-varying point source in a system of two coupled linear diffusion-advection reaction equations: application to surface water pollution[J]. Inverse Problems, 2009, 25(11): 5 009-5 029.

[8]殷鳳蘭,李功勝,賈現正.一個多點源擴散方程的源強識別反問題[J].山東理工大學學報：自然科學版, 2011, 25(2): 01-05.

[9] 谷文娟,李功勝,殷鳳蘭,等.一個時間分數階擴散方程的參數反演問題[J].山東理工大學學報：自然科學版, 2010, 24(6): 22-25.

[10] Li G S, Gu W J, Jia X Z. Numerical inversions for space-dependent diffusion coefficient in the time fractional diffusion equation[J]. Journal of Inverse and Ill-Posed Problems, 2012, 20(3): 339-366.

[11]陸金甫,關治.偏微分方程數值解法[M].北京:清華大學出版社, 2004.

[12] Li G S, Cheng J, Yao D,etal. One dimensional equilibrium model and source parameter determination for soil column experiment[J]. Applied Mathematics and Computation, 2007, 190(2): 1 365-1 374.