婁和忠，李功勝，賈現(xiàn)正

（山東理工大學(xué)理學(xué)院應(yīng)用數(shù)學(xué)所，山東淄博255091）

近年來(lái)，隨著計(jì)算機(jī)技術(shù)的發(fā)展，利用數(shù)值反演確定地下水溶質(zhì)運(yùn)移模型中的彌散系數(shù)或源項(xiàng)系數(shù)等控制參數(shù)的方法越來(lái)越受到地下水工作者們的重視．蘇超偉［1]提出了最佳攝動(dòng)量算法，并對(duì)反應(yīng)擴(kuò)散方程擴(kuò)散系數(shù)的識(shí)別等參數(shù)反演問(wèn)題進(jìn)行了數(shù)值求解．Li等［2－5]利用最佳攝動(dòng)量算法對(duì)一維溶質(zhì)運(yùn)移模型中依賴于空間和時(shí)間變化的源項(xiàng)系數(shù)進(jìn)行了反演研究，并與一維土柱實(shí)驗(yàn)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合，很好地闡釋了實(shí)驗(yàn)結(jié)果．Nie等［6]利用最佳攝動(dòng)量算法對(duì)二維水質(zhì)模型中的綜合擴(kuò)散系數(shù)進(jìn)行了反演，并與檢測(cè)數(shù)據(jù)進(jìn)行了對(duì)比，數(shù)值結(jié)果表明是合理的．閔濤等［7－10]利用有限差分方法對(duì)二維變系數(shù)拋物型方程進(jìn)行了數(shù)值計(jì)算，并且采用遺傳算法，攝動(dòng)量算法以及擬牛頓算法對(duì)模型中的參數(shù)進(jìn)行了反演研究．崔凱等［11]基于同倫映射的思想提出了一種正則參數(shù)依賴于迭代次數(shù)變化的同倫正則化算法，并利用該算法對(duì)一維溶質(zhì)運(yùn)移模型中的多個(gè)參數(shù)進(jìn)行了反演，數(shù)值結(jié)果表明了該算法的有效性．越來(lái)越多的研究表明利用數(shù)值反演的方法確定模型中的未知參數(shù)是非常有效地．總體上看，基于同倫思想對(duì)參數(shù)的反演研究大部分是對(duì)于常數(shù)形式的模型參數(shù)的反演，但當(dāng)未知參數(shù)是隨時(shí)間或空間變化的函數(shù)形式時(shí)，同倫算法的有效性還需進(jìn)一步驗(yàn)證．

1 彌散系數(shù)反演問(wèn)題

本文考慮如下二維對(duì)流彌散方程初邊值問(wèn)題：

初始條件

邊界條件

其中，Ω＝｛（x，y，t）｜0≤x≤1，0≤y≤1，0≤t≤T｝u（x，y，t）為溶質(zhì)濃度，DL（x）＞0為縱向彌散系數(shù)，DT（y）＞0為橫向彌散系數(shù)，v＞0為沿x方向的地下水流速．

如果上述初邊值問(wèn)題中的各個(gè)參數(shù)都已知，利用有限差分方法即可求得其數(shù)值解．然而，在實(shí)際問(wèn)題中，彌散系數(shù)往往是不可知的．選擇觀測(cè)點(diǎn)x＝x0（這里x0是區(qū)域（0，1）內(nèi)的一點(diǎn)），給定某個(gè)終止時(shí)刻t＝T時(shí)的觀測(cè)值作為附加數(shù)據(jù)，即有附加條件

則由（1）－（5）即構(gòu)成了一個(gè)確定彌散系數(shù)DL（x）和DT（y）的參數(shù)反演問(wèn)題．

2 同倫正則化算法

同倫正則化算法屬于顯式迭代算法，通過(guò)使用一種正則參數(shù)依賴于迭代次數(shù)的選取方法，可以避免在迭代過(guò)程中由于試探正則化參數(shù)帶來(lái)的無(wú)效計(jì)算，提高了計(jì)算效率．

本文對(duì)彌散系數(shù)DL和DT的反演采用同倫正則化算法．記D＝（DL，DT），假設(shè)和分別是DL（x）和DT（y）所在空間中的一組基，則DL（x）和DT（y）可以取有限項(xiàng)來(lái)近似如下

和

其中，ai（i＝1，2，…，M）和bj（j＝1，2，…，N）是展開(kāi)項(xiàng)系數(shù)．因此，彌散系數(shù)在已知基函數(shù)的前提下可以寫成D＝（a1，a2，…，aM，b1，b2，…，bN）．

給定彌散系數(shù)的一個(gè)真解D，通過(guò)數(shù)值方法可求得正問(wèn)題的數(shù)值解，記為u（x，y，t；D）．聯(lián)系到附加數(shù)據(jù)θ（y），上述彌散系數(shù)的數(shù)值反演可以化為一個(gè)極小問(wèn)題的求解

其中α＞0為正則參數(shù)．

根據(jù)攝動(dòng)算法的思想，上述極小問(wèn)題的求解又可以轉(zhuǎn)化為對(duì)于給定的初始迭代值DI，求解最佳攝動(dòng)量δDI，即由如下迭代

這里δDI是下述泛函的極小解

基于同倫算法的思想，對(duì)（10）式進(jìn)行不動(dòng)點(diǎn)同倫修正，定義新的目標(biāo)函數(shù)如下

其中，同倫正則參數(shù)λ∈（0，1）采用擬神經(jīng)網(wǎng)絡(luò)函數(shù)定義如下

上式中N0為初值選取參數(shù)，β為下降速率參數(shù)，λ（I）是第I步迭代時(shí)正則參數(shù)的取值．

圖1 同倫正則參數(shù)隨迭代次數(shù)的變化

本文數(shù)值反演中取N0＝1，β＝0．5．圖1給出了同倫正則參數(shù)隨迭代次數(shù)的變化圖像．

注意到上式中的δDI是一個(gè)微小的擾動(dòng)量，將u（x0，y，T；DI＋δDI）在DI處利用多元泰勒公式展開(kāi)得

若對(duì)于（x0，y，T）有Q個(gè)離散點(diǎn)（x0，yq，T）（q＝1，2，…，Q），并且取定‖δDI＝，則得

容易驗(yàn)證［8]，上述泛函極小問(wèn)題的求解等價(jià)于求解以下規(guī)范方程

式中

E為單位矩陣．

其中τ為數(shù)值微分步長(zhǎng)，

下面給出同倫正則化算法的實(shí)施步驟：

步驟1 給定初始迭代向量D0，τ，求解向量η，ξ以及矩陣B；

步驟2 根據(jù)同倫正則化算法求解δDn得到

步驟3 給定收斂精度eps，判斷是否滿足‖δDI‖≤eps．若滿足，則迭代結(jié)束．否則轉(zhuǎn)到式（9）繼續(xù)迭代，直到滿足收斂精度為止．

3 數(shù)值模擬

利用上述介紹的同倫正則化算法對(duì)彌散系數(shù)DL和DT進(jìn)行數(shù)值反演．初邊值問(wèn)題（1）－（4）中取定初始值

邊界條件

附加數(shù)據(jù)取終止時(shí)刻T＝2時(shí)x0＝0．5處的觀測(cè)值．

下面，在式（6）和式（7）中取M＝N＝3，則縱向彌散系數(shù)與橫向彌散系數(shù)可以近似表示為

及

因此，反演參數(shù)可簡(jiǎn)記為

假設(shè)DL（x）＝1＋（y）＝1＋＝1．則反演參數(shù)的真值Dtrue＝（1，0，1，1，0，0．5）．下面應(yīng)用同倫正則化算法對(duì)此彌散系數(shù)進(jìn)行反演重建．

（1）初始迭代值對(duì)反演結(jié)果的影響

取數(shù)值微分步長(zhǎng)τ＝1e－3，迭代收斂精度eps＝1e－6，討論初始迭代值對(duì)算法的影響，計(jì)算結(jié)果見(jiàn)表1，其中Dinv表示反演解，I表示迭代次數(shù)，Err＝‖Dinv－Dtrue表示真解與反演解間的相對(duì)誤差．

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

由表1可知，同倫正則化算法對(duì)初始迭代值的依賴不是很嚴(yán)格．然而，初始迭代值與真解的距離不能太大，否則，算法將失效．

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

取初始迭代值D0＝（0．5，0．5，0．5，0．5，0．5，0．5），迭代收斂精度eps＝le－6，討論微分步長(zhǎng)對(duì)算法的影響，計(jì)算結(jié)果見(jiàn)表2．

由表2可知，數(shù)值微分步長(zhǎng)對(duì)算法的實(shí)現(xiàn)有一定的影響，微分步長(zhǎng)的取值不能太小，不然，算法將不能實(shí)現(xiàn)．

（3）收斂精度對(duì)反演結(jié)果的影響

取初始迭代值D0＝（0．6，0．6，0．6，0．6，0．6，0．6），數(shù)值微分步長(zhǎng)τ＝1e－3，討論收斂精度對(duì)算法的影響，計(jì)算結(jié)果見(jiàn)表3．

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

表3 收斂精度對(duì)反演結(jié)果的影響

由表3可知，收斂精度對(duì)算法的影響很大，收斂精度的取值不能太小．隨著收斂精度的減小，絕對(duì)誤差也逐漸變?。?，收斂精度的取值也不需要太小，否則會(huì)影響反演效率．圖1和圖2分別給出了數(shù)值微分步長(zhǎng)τ＝1e－3，初始值D0＝（0．1，0．1，0．1，0．1，0．1，0．1）時(shí)，縱向彌散系數(shù)DL（x）和橫向彌散系數(shù)DT（y）精確解與反演解得圖像比較．

圖2 DL（x）＝1＋x2時(shí)精確解與反演解的比較

4 結(jié)束語(yǔ)

圖3 DT（y）＝時(shí)精確解與反演解的比較

本文基于同倫正則化算法對(duì)二維對(duì)流－彌散方程中隨時(shí)間變化的縱向彌散系數(shù)和橫向彌散系數(shù)進(jìn)行了反演計(jì)算，研究表明了該算法對(duì)于此類反問(wèn)題參數(shù)反演的有效性．?dāng)?shù)值結(jié)果表明，數(shù)值微分步長(zhǎng)對(duì)反演算法的影響較??；初始迭代值對(duì)算法的影響較大，當(dāng)初始迭代值與精確解的距離較大時(shí)，算法容易失效；收斂精度的取值對(duì)算法的實(shí)現(xiàn)主要是在于對(duì)誤差的影響，收斂精度取的太小反而會(huì)降低計(jì)算效率．

［1] 蘇超偉．偏微分方程逆問(wèn)題的數(shù)值方法及應(yīng)用［M] ．西安：西北工業(yè)大學(xué)出版社，1995

［2] Li G S，Cheng J，Yao D，et al．One－dimensional equilibrium model and source parameter determination for soil－column experiment［J] ．Applied Mathematics and Computation 2007，190（2）：1 365－1 374．

［3] 李功勝，姚德，馬昱，等．一維溶質(zhì)運(yùn)移源（匯）項(xiàng)系數(shù)反演的迭代正則化算法［J] ．地球物理學(xué)報(bào)，2008，51（2）：582－588．

［4] Li G S，Liu J Q，F(xiàn)an X P，et al．A new gradient regularization algorithm for source term inversion in 1Dsolute transportation with final observations［J] ．Applied Mathematics and Computation 2008，196（2）：646－660．

［5] Li G S，Yao D，Jiang H Yet al．Numerical inversion of a time－dependent reaction coefficient in a soil－column infiltrating experiment［J] ．CMES 2011，74（2）：83－107．

［6] Nie N T，Tao J H．Inversion of dispersion coefficient in water quality model using optimal perturbation algorithm［J] ．Applied Mathematics and Mechanics 2009，30（6）：703－712．

［7] 閔濤，張海燕，周宏宇．二維變系數(shù)熱傳導(dǎo)方程初邊值問(wèn)題的交替方向隱格式［J] ．西安工業(yè)大學(xué)學(xué)報(bào)，2007，27（2）：199－204．

［8] 閔濤，張世梅，鄒學(xué)文．二維拋物型方程參數(shù)反演的遺傳算法［J] ．?dāng)?shù)學(xué)雜志，2007，27（3）：348－352．

［9] 閔濤，盧宏鵬，楊曉莉，等．二維變系數(shù)拋物型方程參數(shù)反演的攝動(dòng)量算法［J] ．科技通訊，2010，26（2）：282－287．

［10] 閔濤，盧宏鵬，武苗．二維拋物型方程參數(shù)反演模型的擬牛頓法［J] ．計(jì)算機(jī)工程與應(yīng)用，2010，46（18）：40－42．

［11] 崔凱，李興斯，李寶元．求解非線性反問(wèn)題的大范圍收斂梯度正則化算法［J] ．計(jì)算力學(xué)學(xué)報(bào)，2005，22（4）：415－419．