應用同倫正則化算法反演二維溶質運移模型中的彌散系數

婁和忠，李功勝，賈現正

（山東理工大學理學院應用數學所，山東淄博255091）

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

1 彌散系數反演問題

本文考慮如下二維對流彌散方程初邊值問題：

初始條件

邊界條件

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

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

則由（1）－（5）即構成了一個確定彌散系數DL（x）和DT（y）的參數反演問題．

2 同倫正則化算法

同倫正則化算法屬于顯式迭代算法，通過使用一種正則參數依賴于迭代次數的選取方法，可以避免在迭代過程中由于試探正則化參數帶來的無效計算，提高了計算效率．

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

和

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

給定彌散系數的一個真解D，通過數值方法可求得正問題的數值解，記為u（x，y，t；D）．聯系到附加數據θ（y），上述彌散系數的數值反演可以化為一個極小問題的求解

其中α＞0為正則參數．

根據攝動算法的思想，上述極小問題的求解又可以轉化為對于給定的初始迭代值DI，求解最佳攝動量δDI，即由如下迭代

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

基于同倫算法的思想，對（10）式進行不動點同倫修正，定義新的目標函數如下

其中，同倫正則參數λ∈（0，1）采用擬神經網絡函數定義如下

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

圖1 同倫正則參數隨迭代次數的變化

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

注意到上式中的δDI是一個微小的擾動量，將u（x0，y，T；DI＋δDI）在DI處利用多元泰勒公式展開得

若對于（x0，y，T）有Q個離散點（x0，yq，T）（q＝1，2，…，Q），并且取定‖δDI＝，則得

容易驗證［8]，上述泛函極小問題的求解等價于求解以下規范方程

式中

E為單位矩陣．

其中τ為數值微分步長，

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

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

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

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

3 數值模擬

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

邊界條件

附加數據取終止時刻T＝2時x0＝0．5處的觀測值．

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

及

因此，反演參數可簡記為

假設DL（x）＝1＋（y）＝1＋＝1．則反演參數的真值Dtrue＝（1，0，1，1，0，0．5）．下面應用同倫正則化算法對此彌散系數進行反演重建．

（1）初始迭代值對反演結果的影響

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

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

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

（2）數值微分步長對反演結果的影響

取初始迭代值D0＝（0．5，0．5，0．5，0．5，0．5，0．5），迭代收斂精度eps＝le－6，討論微分步長對算法的影響，計算結果見表2．

由表2可知，數值微分步長對算法的實現有一定的影響，微分步長的取值不能太小，不然，算法將不能實現．

（3）收斂精度對反演結果的影響

取初始迭代值D0＝（0．6，0．6，0．6，0．6，0．6，0．6），數值微分步長τ＝1e－3，討論收斂精度對算法的影響，計算結果見表3．

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

表3 收斂精度對反演結果的影響

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

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

4 結束語

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

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

［1] 蘇超偉．偏微分方程逆問題的數值方法及應用［M] ．西安：西北工業大學出版社，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] 李功勝，姚德，馬昱，等．一維溶質運移源（匯）項系數反演的迭代正則化算法［J] ．地球物理學報，2008，51（2）：582－588．

［4] Li G S，Liu J Q，Fan 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] 閔濤，張海燕，周宏宇．二維變系數熱傳導方程初邊值問題的交替方向隱格式［J] ．西安工業大學學報，2007，27（2）：199－204．

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

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

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

［11] 崔凱，李興斯，李寶元．求解非線性反問題的大范圍收斂梯度正則化算法［J] ．計算力學學報，2005，22（4）：415－419．