楊 振，張耀明, 周愛華，潘月君

(山東理工大學(xué) 理學(xué)院，山東 淄博 255049)

基本解法求解反問題的正則化方法

楊 振，張耀明, 周愛華，潘月君

(山東理工大學(xué) 理學(xué)院，山東 淄博 255049)

針對基本解法在求解反問題時(shí)的病態(tài)特性，將截?cái)嗥娈愔捣纸?TSVD)、Tikhonov正則化方法應(yīng)用于所得病態(tài)系統(tǒng)方程的求解，采用L曲線法和GCV方法確定其正則參數(shù)，并比較了4種組合方法求解的精確性和穩(wěn)定性.數(shù)值算例表明，TSVD方法、Tikhonov正則化方法結(jié)合L曲線法和GCV法可有效地處理反問題中的病態(tài)特性.

基本解法；反問題；截?cái)嗥娈愔捣?；Tikhonov法；L曲線；GCV

數(shù)學(xué)物理反問題的研究興起于20世紀(jì)50年代前后，其研究的主要對象是與勘測、識(shí)別和設(shè)計(jì)等有關(guān)的應(yīng)用問題.在實(shí)際工程應(yīng)用中，有時(shí)由于條件所限，僅僅可以獲得部分邊界上的邊界條件，其余邊界的邊界條件均不可測量.從數(shù)學(xué)的角度來看，這屬于反問題的一種.這類問題由于具有不適定及病態(tài)特性，故其數(shù)值分析具有相當(dāng)大的挑戰(zhàn)性[1-3].邊界元法(BEM)可有效地求解此類問題，但涉及奇異和幾乎奇異積分的處理[4-6]，計(jì)算過程相當(dāng)?shù)胤爆嵑秃臅r(shí).基本解法(MFS)無需對區(qū)域及其邊界進(jìn)行網(wǎng)格劃分，因此無需計(jì)算單元積分，具有計(jì)算精度高、收斂速度快、程序?qū)崿F(xiàn)簡單等優(yōu)良特性，如今已被廣泛應(yīng)用于各種問題的求解，但都不可避免地涉及到病態(tài)線性方程組的處理[6-11].

為了克服基本解法在求解反問題線性系統(tǒng)時(shí)的病態(tài)性，許多學(xué)者提出和發(fā)展了各種方法，其中最具影響的是正則化方法.Chen C S[8]分別用TSVD方法和高斯消去法來求解基本解法中的病態(tài)系統(tǒng)，認(rèn)為TSVD法比高斯消去法優(yōu)越.此后， Feng G R等[9]進(jìn)一步考慮了一般的邊界條件和邊界幾何形狀的問題，得出TSVD法在處理病態(tài)問題時(shí)具有計(jì)算穩(wěn)定、精度高的特點(diǎn).而其他的正則化方法如Tikhonov正則化方法、Landweber迭代法和共軛梯度法等[1-2,11]也有一些應(yīng)用.考慮到目前仍沒有一種適合所有病態(tài)問題解算的最優(yōu)正則化方法，因此對不同的正則化方法進(jìn)行比較研究是一項(xiàng)具有重要意義的工作.

當(dāng)我們選用一種直接正則化方法后，正則化參數(shù)的確定非常重要，它是改善病態(tài)問題解的精確性的關(guān)鍵.目前，正則化參數(shù)的選取準(zhǔn)則有L曲線法、廣義交叉檢驗(yàn)(GCV)準(zhǔn)則等[1-2,11].本文以二維位勢反問題為例，將TSVD、Tikhonov兩種直接正則化方法與L曲線和GCV兩種正則參數(shù)選取方法組合而成的4種不同方法，應(yīng)用于基本解法求解未知邊界條件信息.通過數(shù)值算例，分析比較不同組合方法的計(jì)算效率和精度，以期為實(shí)際應(yīng)用提供合理的選擇依據(jù).

1 二維位勢反問題的基本解法

本文假定Ω是R2中的一個(gè)有界區(qū)域，Γ=?Ω是其邊界.n=(n1,n2)是區(qū)域Ω的邊界Γ在x點(diǎn)處的單位外法向量.

1.1 二維位勢邊界條件識(shí)別反問題

二維位勢問題的控制微分方程為

邊界條件為

式中：u為勢函數(shù)；n為邊界外法向量；Γu, Γq分別是已知u和?u/?n的邊界.

二維位勢問題控制方程的基本解為

1.2 二維位勢反問題的基本解法

基本解法(MFS)是一種無網(wǎng)格法.MFS的基本思想是在邊界之外一定距離處放置一定數(shù)量的虛擬源點(diǎn)，并設(shè)想在每個(gè)源點(diǎn)處存在一個(gè)虛擬的密度函數(shù)值，其在研究區(qū)域上某點(diǎn)產(chǎn)生的影響能夠表示成密度函數(shù)與原問題基本解的線性組合.由于源點(diǎn)分布在物理邊界之外的虛擬邊界上，從而避免了基本解的源點(diǎn)奇異性.目前該方法已被應(yīng)用于各種領(lǐng)域，并取得了很好的效果.

假設(shè)在虛擬邊界?！渖戏植糿個(gè)源點(diǎn)yj,j=1,2,…,n，若在邊界Γ1上選取m個(gè)配置點(diǎn)xi,i=(1,2,…,m).則邊界配置點(diǎn)處的位勢和法向梯度可表示為：

(1)

(2)

這里Φ=(Φ1,Φ2,…,Φn)表示源點(diǎn)的密度函數(shù)值向量，u*(xi,yj)為問題的基本解.

根據(jù)式(1)、式(2)以及已知的邊界條件可得如下一個(gè)線性方程組

(3)

由方程(3)可以計(jì)算出Φ=(Φ1,Φ2,…,Φn)，則由此可通過以下方程計(jì)算出未知邊界上的位勢和通量.

(4)

(5)

一般地，在求解邊界條件識(shí)別反問題時(shí)，方程組(3)往往是不適定的或者高度病態(tài)的，常規(guī)方法如高斯消去法很難求得準(zhǔn)確解.因此在求解方程時(shí)，我們采用正則化方法求解線性方程組(3).

1.3 正則化方法

1.3.1 截?cái)嗥娈愔捣纸?TSVD)

假設(shè)所要求解的方程為

Ax=b

(6)

其中A∈RM×N(M>N),x∈RN,b∈RM，則矩陣A的奇異值分解形式為

(7)

這個(gè)解法用K階矩陣Ak逼近N階矩陣A.下面給出它的形式

(8)

其中K是TSVD方法的截?cái)囗?xiàng)數(shù)，Σk=diag(σ1,σ2,…,σK,0,…,0)，通過這種方式，方程(6)中的矩陣A被AK代替，得到下列方程：

AKx=b

(9)

方程(9)的TSVD解可表示為xK，且

(10)

這里的濾波因子為

1.3.2Tikhonov正則化方法(TR)

Tikhonov正則化方法是一種非常有效而且普遍使用的求解病態(tài)問題的方法，其求解過程是找到一組適定問題的解，使得這組解與原不適定問題比較接近，然后用這組適定問題的解去逼近原不適定問題的解.Tikhonov正則化方法是把正則化泛函

Jα(x)=‖Ax-b‖2+α2‖x‖2,α>0

(11)

的極小元xα作為方程(6)的正則化解.表示成如下形式：

(12)

1.4 正則化參數(shù)的確定

1.4.1L曲線準(zhǔn)則

L曲線準(zhǔn)則的主要思想是通過log-log尺度以反映‖x‖2與‖Ax-b‖2的曲線對比，然后根據(jù)對比得出的結(jié)果決定正則化參數(shù).因?yàn)榍€形狀和“L”的形狀通常比較接近， 所以它被稱為L曲線準(zhǔn)則.

在使用L曲線準(zhǔn)則確定正則化參數(shù)時(shí)，若運(yùn)用Tikhonov正則化方法，則L曲線是由α>0的所有點(diǎn)(‖Axα-b‖,‖xα‖)構(gòu)成的光滑曲線

(13)

若運(yùn)用TSVD法，則L曲線是由一系列點(diǎn)(‖AxK-b‖,‖xK‖)組成的插值樣條曲線.

(14)

對于確定的α>0(或者K=1,2,…,N)，我們解方程(10)(或者(12))就會(huì)得到解xα(或者xK)，以余量?！珹xα-b‖(或者‖AxK-b‖)為橫坐標(biāo)，解的模‖xα‖(或者‖xK‖)為縱坐標(biāo)，可畫出L曲線.找出拐點(diǎn)處的α或者K，即為要找的參數(shù)值.

1.4.2 廣義交叉檢驗(yàn)(GCV)準(zhǔn)則

GCV準(zhǔn)則是用正則化參數(shù)α(或者K)作為參變量，來求解GCV函數(shù)的最小值，只要得到GCV函數(shù)的極小值，那么對應(yīng)α(或者K)就是我們要確定的最優(yōu)正則化參數(shù).其計(jì)算公式為

(15)

或者

(16)

2 數(shù)值算例

考慮兩個(gè)二維位勢反問題的數(shù)值算例，采用4種不同的正則化組合方法對算例進(jìn)行求解.為了對比不同方法數(shù)值解的準(zhǔn)確性，定義如下平均相對誤差

(17)

例1 如圖1所示方形閉域熱流問題，邊長為2.精確解為u=y2-x2，邊界條件為：邊AB、BC上的溫度和熱流已知，邊CD、DA上所有邊界條件未知.

如圖1所示，在距離真實(shí)邊界距離為d的地方取一虛擬邊界?！洌蛊渑c熱流區(qū)域邊界相似.在虛擬邊界四邊各分布10個(gè)源點(diǎn)，共40個(gè)源點(diǎn)，如圖2所示.為了考察正則化方法對于位勢反問題求解的有效性，我們分別用4種組合方法對本例進(jìn)行求解.我們?nèi)≡袋c(diǎn)分布在虛擬邊界Γ′上，選取不同的虛實(shí)邊界距離d分別為5和10.圖3～圖 6給出了d取不同值時(shí)，未知邊界上各點(diǎn)溫度和熱流的相對誤差，可以看出，本文給出的4種組合方法在未知邊界各點(diǎn)求得的數(shù)值解都獲得了較高精度.

圖1 方形閉域熱流 圖2 節(jié)點(diǎn)編號(hào)

圖3 未知邊界的溫度誤差(d=5)

圖4 未知邊界的熱流誤差(d=5)

圖5 未知邊界的溫度誤差(d=10)

圖6 未知邊界的熱流誤差(d=10)

例2 如圖7所示圓形閉域熱流問題，半徑為1.精確解為u=x2-y2，我們考慮的邊界條件為：1/4邊界上的溫度和熱流已知，另外3/4邊界上的邊界條件未知.此時(shí)我們所得方程是不適定的，因此常規(guī)方法很難求解.

如圖7所示，在距離真實(shí)邊界距離為d的地方取一虛擬邊界Γ′，使其與熱流區(qū)域邊界相似.并在虛擬邊界上均勻地布置100個(gè)源點(diǎn).同樣我們分別用4種組合方法對本例進(jìn)行求解.選取不同的虛實(shí)邊界距離d分別為5和10.表1、表2給出了d取不同值時(shí)，未知邊界上溫度和熱流的平均相對誤差.表1和表2更準(zhǔn)確地反映出本文給出的4種組合方法計(jì)算結(jié)果差異，也不難看出，MFS結(jié)合4種組合方法都能夠有效地求解不同邊界幾何形狀的位勢邊界條件反問題.

圖7 圓形閉域熱流

表1 四種方法參數(shù)選取及平均相對誤差比較結(jié)果(d=5)

正則化方法參數(shù)(α/K)平均相對誤差/%溫度熱流TSVD-LC173．400991×10-61．709062×10-6TSVD-GCV122．492231×10-53．790088×10-5TR-LC2．44×10-113．268636×10-61．602206×10-6TR-GCV2．01×10-68．815173×10-51．121432×10-4

表2 四種方法參數(shù)選取及平均相對誤差比較結(jié)果(d=10)

正則化方法參數(shù)(α/K)平均相對誤差/%溫度熱流TSVD-LC131．946931×10-79．724293×10-8TSVD-GCV122．417499×10-73．138218×10-7TR-LC3．37×10-111．577152×10-76．523213×10-8TR-GCV3．91×10-111．866861×10-79．110862×10-8

3 結(jié)束語

本文針對基本解法求解反問題時(shí)的病態(tài)特性，應(yīng)用TSVD方法和Tikhonov正則化方法處理所得病態(tài)系統(tǒng)方程，L曲線法和GCV法確定其正則化參數(shù)，并比較了不同正則化方法的計(jì)算精度.數(shù)值算例表明，針對不同邊界條件以及邊界幾何形狀，MFS結(jié)合TSVD方法和Tikhonov正則化方法都能夠有效求解二維位勢反問題，即使已知邊界信息非常有限，依然能夠獲得穩(wěn)定、高精度的數(shù)值解.

[1]肖庭延. 反問題的數(shù)值解法[M]. 北京：科學(xué)出版社, 2003.

[2]劉繼軍. 不適定問題的正則化方法及應(yīng)用[M]. 北京：科學(xué)出版社, 2005.

[3]HansenPC.Rank-deficientanddiscreteill-poseedproblems:numericalaspectsoflinearinversion[M].Philadelphia:SIAM,1998.

[4]張耀明,谷巖,陳正宗. 位勢邊界元法中的邊界層效應(yīng)與薄體結(jié)構(gòu)[J]. 力學(xué)學(xué)報(bào)，2010，42(2):219-227.

[5]ZhangYM,SunCL.Ageneralalgorithmforthenumericalevaluationofnearlysingularboundaryintegralsintheequivalentnon-singularBIEswithindirectunknowns[J].JournaloftheChineseInstituteofEngineers, 2008, 31(3):437-447.

[6]孫煥純. 無奇異邊界元法[M]. 大連：大連理工大學(xué)出版社, 1999.

[7]KupradzeVD，AleksidzeMA.Themethodoffunctionalequationsfortheapproximatesolutionofcertainboundaryvalueproblems[J].USSRComput.MathPhy，1964,4:82-126.

[8]ChenCS,HokwonA.Somecommentsontheill-conditioningofthemethodoffundamentalsolutions[J].EngineeringAnalysiswithBoundaryElements,2006,30:405-410

[9]FengGR,LiM,ChenCS.Ontheill-conditioningoftheMFSforirregularboundarydatawithsufficientregularity[J].EngineeringAnalysiswithBoundaryElements,2014,41:98-102.

[10]TikhonovAN,ArseninVY.SolutionsofILL-posedproblem[M].NewYork:JohnWileyandSons，1977.

[11]WeiT,HonYC,LingLV.MethodoffundamentalsolutionswithregularizationtechniquesforCauchyproblemsofellipticoperators[J].EngineeringAnalysiswithBoundaryElements,2007,31:373-385.

(編輯：郝秀清)

The regularization method of the MFS for the inverse problem

YANG Zhen, ZHANG Yao-ming， ZHOU Ai-hua， PAN Yue-jun

(School of Science, Shandong University of Technology, Zibo 255049, China)

In order to resolve ill-conditioned problems existing in the regularization MFS for the 2D boundary conditions identification potential problems, suitable regularization methods are needed. The truncated singular value decomposition(TSVD)method and Tikhonov regularization method are used to solve linear systems with a large number of conditions respectively. The L-curve and generalized cross validation (GCV) methods are employed to determine the optimal regularization parameters. Furthermore, the accuracy and robustness of regularization solution for four combined methods are investigated． Numerical results show that TSVD method and Tikhonov method can effectively solve the ill posed system caused by inverse problems. Through applying to L-curve method and GCV method, continuous regularization parameter for Tikhonov method can be confirmed reasonably.

MFS; inverse problems; truncated singular value decomposition; Tikhonov method;L-curve; GCV

2015-01-11

山東省自然科學(xué)基金重點(diǎn)項(xiàng)目(ZR2010AZ003)

楊振, 男, librayz@126.com; 通信作者：張耀明，男，zymfc@163.com

1672-6197(2015)06-0020-05

O342

A