時(shí)間分?jǐn)?shù)階二維對(duì)流擴(kuò)散方程多點(diǎn)源強(qiáng)的數(shù)值反演

李慧玲， 李功勝， 賈現(xiàn)正， 池光勝

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

分布式參數(shù)系統(tǒng)模型的建立對(duì)精確描述和有效控制許多物理和工程問題起關(guān)鍵作用.源項(xiàng)的估計(jì)或識(shí)別問題在環(huán)境科學(xué)和工業(yè)應(yīng)用等領(lǐng)域顯得尤為重要.譬如，在環(huán)境水力學(xué)領(lǐng)域，如何尋找河流、城市水環(huán)境和湖泊的污染源；在化學(xué)反應(yīng)過程、微波加熱過程和工業(yè)設(shè)計(jì)與制造領(lǐng)域，如何探測(cè)未知的發(fā)熱源等.

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

本文將探討時(shí)間分?jǐn)?shù)階二維對(duì)流擴(kuò)散方程中確定多個(gè)點(diǎn)源的反演問題.對(duì)于這類分?jǐn)?shù)階擴(kuò)散方程點(diǎn)源識(shí)別反問題的研究尚未見有文獻(xiàn)報(bào)道.文中將首先給出正問題及其數(shù)值求解的差分格式，然后引入最佳攝動(dòng)量正則化算法，并對(duì)多個(gè)點(diǎn)源強(qiáng)度進(jìn)行精確數(shù)據(jù)和擾動(dòng)數(shù)據(jù)條件下的數(shù)值反演.反演結(jié)果表明對(duì)于時(shí)間分?jǐn)?shù)階二維對(duì)流擴(kuò)散問題，應(yīng)用最佳攝動(dòng)量正則化算法可以實(shí)現(xiàn)對(duì)多個(gè)點(diǎn)源強(qiáng)度的數(shù)值確定.

1 正問題及其數(shù)值求解

考慮矩形域Ω=(0，l1)×(0,l2)上的時(shí)間分?jǐn)?shù)階二維擴(kuò)散方程

(1)

其中(x,y)∈Ω,t>0;D>0是擴(kuò)散系數(shù)，v>0是沿著X軸方向的對(duì)流速度，g(x,y,t)是線性源項(xiàng)，而α∈(0,1)是時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù).當(dāng)α=1時(shí),方程(1)即為通常的整數(shù)階二維對(duì)流擴(kuò)散方程.對(duì)于帶非零源項(xiàng)的方程(1),給定零初邊值條件:

(2)

這樣,由方程(1)及初邊值條件(2),就構(gòu)成一個(gè)二維時(shí)間分?jǐn)?shù)階對(duì)流擴(kuò)散的正問題. 下面給出這個(gè)正問題求解的一個(gè)有限差分格式.

1.1 正問題求解的差分格式

O(τ).

(3)

對(duì)于方程中的整數(shù)階導(dǎo)數(shù)項(xiàng)，按照通常的差分離散方法(參見文獻(xiàn)[11])，即有

(4)

(5)

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

當(dāng)n=0時(shí)

(6)

當(dāng)n>0時(shí)

(7)

1.2 數(shù)值試驗(yàn)

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

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

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

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

2 多點(diǎn)源項(xiàng)識(shí)別反問題與反演算法

2.1 反問題的提出

很多實(shí)際問題中,方程(1)中的源項(xiàng)是點(diǎn)源分布,且一般具有形式

(8)

其中：q為點(diǎn)源的個(gè)數(shù)；(xs,ys)為點(diǎn)源的位置；Qs為相應(yīng)的點(diǎn)源強(qiáng)度(單位時(shí)間排放量)；δ為狄拉克函數(shù).

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

VT

(9)

這樣，由附加數(shù)據(jù)(9)聯(lián)合正問題(1)-(2)構(gòu)成了一個(gè)確定源強(qiáng)度Qs(s=1,2,…，q)的多點(diǎn)源強(qiáng)識(shí)別反問題.下面給出最佳攝動(dòng)量正則化算法，并對(duì)上述源強(qiáng)識(shí)別反問題進(jìn)行數(shù)值反演模擬.

2.2 最佳攝動(dòng)量正則化算法

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

(10)

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

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

(11)

且δQn是下述目標(biāo)函數(shù)的極小值

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),

略去高階項(xiàng)，則目標(biāo)函數(shù)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)=

λ為數(shù)值微分步長(zhǎng). 于是

μ‖δQn‖2

(14)

再根據(jù)最小二乘法的思想，求解minF(δQn)相當(dāng)于求解規(guī)范方程

[ATA+μI]δQn=ATB

(15)

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

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

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

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

3 數(shù)值反演

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

3.1 兩個(gè)點(diǎn)源的情形

設(shè)附加數(shù)據(jù)在(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個(gè)點(diǎn)處取得.假設(shè)在(0.2, 0.2), (0.8, 0.8)處有兩個(gè)點(diǎn)源，其源強(qiáng)分別為2和10，即源強(qiáng)真值為Q=[2,10].

(1)正則參數(shù)對(duì)反演結(jié)果的影響

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

表4 正則參數(shù)對(duì)反演結(jié)果的影響

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

(2)微分階數(shù)對(duì)反演結(jié)果的影響

取數(shù)值微分步長(zhǎng)λ=0.1，初始迭代值Q0=[0,0]，計(jì)算結(jié)果見表5.

表5 微分階數(shù)對(duì)反演結(jié)果的影響

通過表5可以看出微分階數(shù)對(duì)反演結(jié)果的影響不大.

(3)初始迭代對(duì)反演結(jié)果的影響

取微分階數(shù)α=0.6，數(shù)值微分步長(zhǎng)λ=0.1，計(jì)算結(jié)果見表6.

表6 初始迭代對(duì)反演結(jié)果的影響

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

(4)數(shù)值微分步長(zhǎng)對(duì)反演結(jié)果的影響

取微分階數(shù)α=0.6，初始迭代值Q0=[0,0]，計(jì)算結(jié)果見表7.

表7 數(shù)值微分步長(zhǎng)對(duì)反演結(jié)果的影響

通過表7可以看出數(shù)值微分步長(zhǎng)對(duì)反演結(jié)果也幾乎沒有影響.

3.2 三個(gè)點(diǎn)源的情形

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

(1)擾動(dòng)水平對(duì)反演結(jié)果的影響

仍取微分階數(shù)α=0.6，數(shù)值微分步長(zhǎng)λ=0.1，分別取擾動(dòng)水平δ=10%,5%,1%進(jìn)行反演計(jì)算,10次反演的平均值列于表8.

表8 擾動(dòng)水平對(duì)反演結(jié)果的影響

通過表8可以看出隨著數(shù)據(jù)擾動(dòng)水平的減小，反演結(jié)果精度越來越高.

(2)分?jǐn)?shù)微分階數(shù)對(duì)反演結(jié)果的影響

設(shè)擾動(dòng)水平δ=1%，考察分?jǐn)?shù)微分階數(shù)對(duì)反演算法的影響，計(jì)算結(jié)果見表9.

表9 給定擾動(dòng)水平下分?jǐn)?shù)微分階數(shù)對(duì)反演結(jié)果的影響

由表9可以看出,對(duì)于給定的擾動(dòng)數(shù)據(jù),分?jǐn)?shù)微分階數(shù)對(duì)反演結(jié)果的影響不大.

3.3 四個(gè)點(diǎn)源的情形

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

(1)附加數(shù)據(jù)的選取對(duì)反演結(jié)果的影響

取微分階數(shù)α=0.7，數(shù)值微分步長(zhǎng)λ=0.1，附加數(shù)據(jù)選取x=0.1直線段上各點(diǎn)(0.1,yk)(k=1,…，K)處的觀測(cè)值.應(yīng)用精確數(shù)據(jù)進(jìn)行反演，對(duì)應(yīng)于不同觀測(cè)點(diǎn)的計(jì)算結(jié)果見表10，其中y=[·]表示在x=0.1線上所取觀測(cè)點(diǎn)的縱坐標(biāo).

表10 附加數(shù)據(jù)的選取對(duì)反演結(jié)果的影響

通過表10可以看出，附加數(shù)據(jù)對(duì)反演結(jié)果具有一定的影響.若獲得附加數(shù)據(jù)的觀測(cè)點(diǎn)少于4個(gè),則反演計(jì)算難以得到理想結(jié)果.但當(dāng)選取4個(gè)以上觀測(cè)數(shù)據(jù)點(diǎn)時(shí),反演結(jié)果相差并不大.

(2) 附加數(shù)據(jù)的擾動(dòng)對(duì)反演結(jié)果的影響

仍在x=0.1直線段上選取觀測(cè)點(diǎn),不妨設(shè)有9個(gè)觀測(cè)點(diǎn),即x=0.1,y=[0.2,0.3,…，0.9], 考察附加數(shù)據(jù)的隨機(jī)擾動(dòng)對(duì)反演算法的影響.十次反演計(jì)算的平均結(jié)果列于表11.

表11 附加數(shù)據(jù)擾動(dòng)對(duì)反演結(jié)果的影響

通過表11可以看出,在四個(gè)點(diǎn)源的情形,數(shù)據(jù)擾動(dòng)對(duì)反演算法影響較大,這說明此時(shí)的反演問題病態(tài)性較重.不過,隨著擾動(dòng)水平的減小,反演誤差也變得越來越小.

4 結(jié)束語

本文探討了二維時(shí)間分?jǐn)?shù)階對(duì)流擴(kuò)散方程中多點(diǎn)源強(qiáng)度的識(shí)別反問題.基于正問題求解的差分格式,應(yīng)用最佳攝動(dòng)量算法進(jìn)行了有效的數(shù)值反演.反演結(jié)果表明這類源項(xiàng)反問題的病態(tài)性并不強(qiáng),特別當(dāng)點(diǎn)源個(gè)數(shù)較少時(shí),完全可以采用不帶顯式正則化項(xiàng)的最佳攝動(dòng)量算法進(jìn)行反演確定.對(duì)于四個(gè)點(diǎn)源的情形,只要給定四個(gè)附加數(shù)據(jù)點(diǎn)就可以實(shí)施模擬計(jì)算,但此時(shí)附加數(shù)據(jù)的擾動(dòng)對(duì)反演算法影響較大,這也說明在點(diǎn)源個(gè)數(shù)多于四個(gè)的情形,反演問題的病態(tài)性增大,需要尋找更為有效的反演方法才能獲得更穩(wěn)定的結(jié)果.

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

[2]閔濤,周孝德,張世梅,等.對(duì)流-擴(kuò)散方程源項(xiàng)識(shí)別反問題的遺傳算法[J].水動(dòng)力學(xué)研究與進(jìn)展,2004,19(4):520-524.

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

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

[5]張小明,王澤文,王燕等.一類熱傳導(dǎo)方程多點(diǎn)源反演的唯一性和穩(wěn)定性[J].華東理工大學(xué)學(xué)報(bào)(自然科學(xué)版) ,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]殷鳳蘭,李功勝,賈現(xiàn)正.一個(gè)多點(diǎn)源擴(kuò)散方程的源強(qiáng)識(shí)別反問題[J].山東理工大學(xué)學(xué)報(bào)：自然科學(xué)版, 2011, 25(2): 01-05.

[9] 谷文娟,李功勝,殷鳳蘭,等.一個(gè)時(shí)間分?jǐn)?shù)階擴(kuò)散方程的參數(shù)反演問題[J].山東理工大學(xué)學(xué)報(bào)：自然科學(xué)版, 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]陸金甫,關(guān)治.偏微分方程數(shù)值解法[M].北京:清華大學(xué)出版社, 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.