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

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

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

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

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

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

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

(1)

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

(2)

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

1.1 正問題求解的差分格式

O(τ).

(3)

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

(4)

(5)

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

當(dāng)n=0時

(6)

當(dāng)n>0時

(7)

1.2 數(shù)值試驗

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

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

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

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

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

2.1 反問題的提出

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

(8)

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

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

VT

(9)

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

2.2 最佳攝動量正則化算法

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

(10)

根據(jù)最佳攝動量算法(參見文獻[8-10, 12])，上述極小問題(10)的求解又轉(zhuǎn)化為對于給定的Qn,通過求解最佳攝動量δQn進而確定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),

略去高階項，則目標(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ù)值微分步長. 于是

μ‖δQn‖2

(14)

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

[ATA+μI]δQn=ATB

(15)

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

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

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

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

3 數(shù)值反演

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

3.1 兩個點源的情形

設(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個點處取得.假設(shè)在(0.2, 0.2), (0.8, 0.8)處有兩個點源，其源強分別為2和10，即源強真值為Q=[2,10].

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 三個點源的情形

設(shè)在點(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)擾動水平對反演結(jié)果的影響

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

表8 擾動水平對反演結(jié)果的影響

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

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

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

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

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

3.3 四個點源的情形

設(shè)在(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)附加數(shù)據(jù)的選取對反演結(jié)果的影響

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

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

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

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

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

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

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

4 結(jié)束語

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

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

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

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

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

[5]張小明,王澤文,王燕等.一類熱傳導(dǎo)方程多點源反演的唯一性和穩(wěn)定性[J].華東理工大學(xué)學(xué)報(自然科學(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)正.一個多點源擴散方程的源強識別反問題[J].山東理工大學(xué)學(xué)報：自然科學(xué)版, 2011, 25(2): 01-05.

[9] 谷文娟,李功勝,殷鳳蘭,等.一個時間分?jǐn)?shù)階擴散方程的參數(shù)反演問題[J].山東理工大學(xué)學(xué)報：自然科學(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.