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

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

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

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

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

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

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

(1)

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

(2)

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

1.1 正問題求解的差分格式

O(τ).

(3)

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

(4)

(5)

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

當n=0時

(6)

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

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

2.1 反問題的提出

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

(8)

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

假設已經(jīng)知道點源的個數(shù)及其位置坐標，那么為了確定各個點源的強度值，需要補充關于污染物濃度分布的附加條件，并聯(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)構成了一個確定源強度Qs(s=1,2,…，q)的多點源強識別反問題.下面給出最佳攝動量正則化算法，并對上述源強識別反問題進行數(shù)值反演模擬.

2.2 最佳攝動量正則化算法

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

(10)

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

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

(11)

且δQn是下述目標函數(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),

略去高階項，則目標函數(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)相當于求解規(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，再轉到步驟1繼續(xù)進行.

3 數(shù)值反演

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

3.1 兩個點源的情形

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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)擾動水平對反演結果的影響

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結束語

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

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

[9] 谷文娟,李功勝,殷鳳蘭,等.一個時間分數(shù)階擴散方程的參數(shù)反演問題[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]陸金甫,關治.偏微分方程數(shù)值解法[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.