胡宇清

(東南大學數(shù)學系，211100，南京)

0 引言

在科學研究中經(jīng)常要通過間接觀測來探索位于不可達、不可觸之處的物質(zhì)的變化規(guī)律;生產(chǎn)中經(jīng)常要根據(jù)特定的功能對產(chǎn)品進行設(shè)計，或按照某種目的對流程進行控制。簡而言之，就是由果索因，而這就是所謂的反問題。反問題大量出現(xiàn)于醫(yī)學成像、地質(zhì)探測、非損傷性探測、雷達遙感等領(lǐng)域。反問題不僅是非線性的，而且從數(shù)值計算的角度看是不適定的。

熱成像、靜電成像等很多物理問題在數(shù)學上被看做Laplace方程的逆邊值問題，本文考慮這類問題的廣義阻尼邊值條件(GIBC)。目前對廣義阻尼邊界條件問題已有很多研究[1－3]。對于雙連通區(qū)域已知外部邊界上一些柯西數(shù)據(jù)重構(gòu)內(nèi)部邊界形狀的反問題通常采用文獻[4]中提出的邊界積分方程方法。邊界積分方程方法已被廣泛應用，文獻[5，6]中作者應用該方法重構(gòu)了腐蝕探測問題中被腐蝕邊界的形狀以及阻尼系數(shù);文獻[7]中作者對已知的腔體內(nèi)部一條曲線上若干散射場的測量數(shù)據(jù)，運用該方法重構(gòu)了外部腔體的形狀。本文基于文獻[4]中提出的方法通過格林公式定理推導出邊界積分方程組，從而將要解決的重構(gòu)邊界形狀反問題轉(zhuǎn)化為了求解非線性積分方程組的問題，再用正則化的迭代方法最終重構(gòu)出邊界形狀。

1 問題的敘述

假設(shè)雙連通有界區(qū)域D?R2，其邊界由2個非連接的光滑若當閉曲線Γ0和Γ1組成，使得?D∶ =Γ0∪Γ1且Γ0∩Γ1=φ，其中Γ1是環(huán)形區(qū)域D的內(nèi)邊界。定義υ為區(qū)域D的邊界單位外法向量?？紤]如下邊值問題:對給定的函數(shù)(Γ0)，u∈H2(D)滿足

其中:divΓ1和 gradΓ1分別為表面散度和梯度，λ∈C1(Γ1)非負且不恒為0，μ∈C1(Γ1)為正。

假定邊界Γ0上Dirichlet數(shù)據(jù)f對應的Neumann數(shù)據(jù)

是已知的。要考慮的反問題為:假設(shè)已知邊界Γ1上阻尼系數(shù)，由Γ0上的一組柯西數(shù)據(jù)

1(Γ0)×H2(Γ0)重構(gòu)邊界Γ1的形狀。在二維情況下，邊界條件中的Laplace-Beltrami微分算子有表達式為切向?qū)?shù)，s為弧長。盡管這類反問題解的唯一性是未知的，但本文提出的數(shù)值反演方法仍然能很好地得到重建結(jié)果。

2 非線性積分方程

先從格林公式出發(fā)得到與原來反問題等價的非線性積分方程。根據(jù)Laplace方程的基本解

它們限制在部分邊界上的形式為

其中，j，k=0，1。

對邊值問題(1)的解u∈H2(D)，令α∶ =u|Γ1并簡記

由格林公式[8]得到

在上式中令x由D的內(nèi)部分別接近邊界Γ0和Γ1，根據(jù)雙層位勢在邊界上的跳躍關(guān)系可得積分方程

和

因此，要求的邊界形狀重構(gòu)問題就轉(zhuǎn)化為了由上式2個非線性方程組求解Γ1和α的問題。式(5)、式(6)是典型的非線性不適定積分方程組，這2個式子對于邊界Γ1都是非線性的，對于α都是線性的。求解式(5)和式(6)的步驟是首先利用Frechet導數(shù)把非線性方程組近似表示為其線性形式，再用Tikhonov正則化方法處理問題的不適定性。

有3種迭代法求解非線性積分方程組(5)、(6)。

1)從式(5)中由線性方程解出α，代入式(3)得到β，然后關(guān)于未知邊界Γ1線性化方程(6)，從而迭代更新得近似邊界形狀。

2)從式(6)出發(fā)，計算出α再代入式(3)得到β，對邊界Γ1線性化方程(5)迭代更新得到近似邊界形狀。

3)將式(5)、式(6)關(guān)于邊界Γ1和α同時線性化，解這2個方程，得到關(guān)于Γ1和α的更新迭代。文中將采用第3種迭代方法重構(gòu)得到近似邊界形狀。

3 積分方程離散和迭代法

不失一般性，假定C2邊界Γ0，Γ1的參數(shù)化形式分別為

其中，z0∶R→R2，z1∶R→R2為以2π 為周期的 C2光滑函數(shù)。對任意一個向量a=(a1，a2)定義a⊥=(a2，－a1)為向量a順時針旋轉(zhuǎn)90°。令ψ=φ?zj，j=0，1，根據(jù)單雙層位勢算子定義得對應的離散化的積分算子形式分別為

式(9)中第2項的對角線元素為

將式(5)、式(6)關(guān)于未知數(shù)Γ1和α完全線性化得到

1)給出未知邊界Γ1的初始參數(shù)化估計z1。已知邊界Γ1上阻尼系數(shù)λ，μ以及Γ0上的柯西數(shù)據(jù)f，g利用式(5)得到α，從而由式(3)計算出β。

2)將α，β以及邊界初始估計z1代入式(10)和式(11)，解線性方程組得γ，ζ，從而得到更新α+ γ，z1+ ζ。

3)重復第2步直到滿足迭代終止條件。下面給出Frechet導數(shù)的計算形式。

式(14)第1項和第3項中核的對角線元素為

第2項和第4項核的對角線元素可由式(8)得到。

ω關(guān)于z1的Frechet導數(shù)為

式(12)～式(16)中 t∈[0，2π]。下一節(jié)就將上述迭代法給出具體算例。

4 數(shù)值算例

根據(jù)文獻[3]中介紹的利用單層位勢在邊界上的跳躍關(guān)系由已知邊界上的Dirichlet數(shù)據(jù)算出Neumann數(shù)據(jù)，以下重構(gòu)邊界形狀的數(shù)值算例都只用一組柯西數(shù)據(jù)，然后用積分方程方法重構(gòu)內(nèi)部邊界形狀。

將邊界Γ1的更新記做

其中，q為非負函數(shù)。假設(shè)q可由如下三角多項式的形式來近似

在所有數(shù)值算例中外部的已知曲線Γ0為圓心在原點半徑為0.9的圓，Γ1的初始估計為圓心在原點半徑為0.8的圓，配置點個數(shù)為n=32，Γ1上三角多項式的展開系數(shù)取m=6。記αγ，αζ分別為γ以及邊界更新ζ的正則化參數(shù)，it記為迭代步數(shù)。為了檢驗該迭代算法的穩(wěn)定性，給測量數(shù)據(jù)g加上誤差，誤差數(shù)據(jù)由如下式子給出

其中，η是正態(tài)分布隨機變量，δ是相對噪音水平。

例1:考慮橢圓邊界，其參數(shù)化形式為

給定Dirichlet數(shù)據(jù)f(z0(t))=1。運用上述的正則化迭代法得到圖1的重構(gòu)結(jié)果。

圖1中左圖選取的正則化參數(shù)為αγ=10－2，αζ=10－1，迭代步數(shù)it=7。圖1中右圖選取的正則化參數(shù)為 αγ=10－2，αζ=10－2，迭代步數(shù) it=9。圖1中邊界Γ1上阻尼系數(shù)取常數(shù)λ=2，μ=2，下面考慮阻尼系數(shù)為連續(xù)函數(shù)的情況。取λ=，Dirichlet數(shù)據(jù)仍為 f(z0(t))=1，得到圖2所示重構(gòu)結(jié)果。

圖1 用精確測量數(shù)據(jù)和6%擾動數(shù)據(jù)重構(gòu)橢圓邊界形狀

圖2 用精確測量數(shù)據(jù)和6%擾動數(shù)據(jù)重構(gòu)橢圓邊界形狀

圖2中左圖選取的正則化參數(shù)為αγ=10－2，αζ=10－1，迭代步數(shù)it=9。圖2中右圖選取的正則化參數(shù)為 αγ=10－2，αζ=10－2，迭代步數(shù) it=7。

例2:考慮梨邊界，其參數(shù)化形式為

給定Dirichlet數(shù)據(jù) f(z0(t))=1，Γ1上的阻尼系數(shù)為，同樣用精確數(shù)據(jù)和擾動數(shù)據(jù)可以得到圖3所示的結(jié)果。

圖3 用精確測量數(shù)據(jù)和3%擾動數(shù)據(jù)重構(gòu)梨邊界形狀

圖3中2個圖選取的正則化參數(shù)均為αγ=9×10－1，αζ=10－2，迭代步數(shù) it=6?？紤]對不同的Dirichlet數(shù)據(jù)f(z0(t))=1+sin2t重構(gòu)梨邊界形狀，Γ1上的阻尼系數(shù)仍為=1+cos4t，得到圖4所示的重構(gòu)結(jié)果。

圖4中2個圖選取的正則化參數(shù)均為αγ=10－1，αζ=5 ×10－1，迭代步數(shù) it=6。這個算例說明對邊界Γ0上不同的柯西數(shù)據(jù)提出的正則化迭代算法仍然是有效的。

例3:考慮蘋果邊界，其參數(shù)化形式為

取定 Dirichlet數(shù)據(jù) f(z0(t))=1+sin2t，Γ1上的阻尼系數(shù)為，得到圖 5所示的重構(gòu)結(jié)果。

圖4 用精確測量數(shù)據(jù)和3%擾動數(shù)據(jù)重構(gòu)梨邊界形狀

圖5 用精確測量數(shù)據(jù)和3%擾動數(shù)據(jù)重構(gòu)蘋果邊界形狀

圖5中左圖選取的正則化參數(shù)為αγ=10－3，αζ=6 ×10－2，迭代步數(shù) it=13。圖 5 中右圖選取的正則化參數(shù)為 αγ=10－3，αζ=10－3，迭代步數(shù) it=15。

通過上述3個數(shù)值例子可以看出，本文提出的重構(gòu)內(nèi)部邊界形狀的非線性積分方程法是有效的，無論Γ1上取不同阻尼系數(shù)，還是Γ0上取不同的Dirichlet數(shù)據(jù)，該算法都能很好的重構(gòu)內(nèi)部邊界形狀，并且對數(shù)據(jù)誤差是穩(wěn)定的。

致謝:感謝導師劉繼軍教授在本文完成中的討論和提出的有效建議。

[1] Haddar H，Joly P，Nguyen H M.Generalized impedance boundary conditions for scattering problems from strongly absorbing obstacles:the case of maxwell's equations[J].Math Models and Methods in Appl Sci，2008，18:1787－127.

[2] Bourgeois L，Chaulet N，Haddar H.Stable reconstruction of generalized impedance boundary conditions[J].Inverse Problems，2011，27:095002 －26.

[3] Cakoni F，Kress R.Integral equation methods for the inverse obstacle problem with generalized impedance boundary condition[J].Inverse Problems，2013，29:015005－19.

[4] Kress R，Rundell W.Nonlinear integral equations and the iterative solution for an inverse boundary value problem[J].Inverse Problems，2005，21:1207 －1223.

[5] Cakoni F，Kress R，Schuft C.Integral equations for shape and impedance reconstruction in cottosion detection[J].Inverse Problems，2010，26:095012 －24.

[6] Cakoni F，Kress R，Schuft C.Simultaneous reconstruction of shape and impedance in corrosion detection[J].Methods Appl Anal，2010，17:331 －462.

[7] Qin H H，Cakoni F.Nonlinear integral equations for shape reconstruction in the inverse interior scattering problem[J].Inverse Problems，2011，27:035005 － 17.

[8] Colton D，Kress R.Inverse Acoustic and ElEctromagnetic Scattering Theory[M].Berlin:Springer，1998.