【作者】 张磊；

【导师】 马富明；

【作者基本信息】 吉林大学 ， 计算数学， 2013， 博士

【摘要】 无界粗糙表面散射与反散射问题在科学研究和工程实际中,特别在地学遥感,目标识别,光学衍射等领域中有十分重要的应用.声波和电磁波的散射与反散射问题是现代科学领域中研究的热点问题之一.散射问题是指对于给定的入射场及散射体,确定波场的散射特性,如散射场或散射场在无穷远处的性态；而反散射问题是根据给定的入射场和散射场(或远场)的信息,求出散射体的边界或者相关的物理参数.求解散射问题常用的方法有积分方程方法,有限元法和无限元法等数值方法.求解反散射问题的主要数值方法有迭代法,优化方法和抽样探测法等.本文主要针对锥形波入射情形的无界粗糙表面散射和反散射问题进行研究.首先,从模型问题出发,严格推导了锥形波入射情形的无界粗糙表面散射问题的相应积分方程.同时证明了积分算子的一些相关性质.然后,根据锥形波入射情形下相应积分方程的具体特点,提出了快速多极法与正则化共轭梯度法相结合的数值方法求解锥形波入射情形的无界粗糙表面散射问题.最后,讨论了锥形入射波情形下完全导体(声软)无界粗糙表面的形状重构问题.具体工作如下：Ⅰ.锥形波入射无界粗糙表面散射问题边界积分方程方法设Γ：={(x,f(x))∈R2|x∈R}表示无界粗糙表面,Ω：={r=(x,z)∈R2|z>f(x)}是Γ上部传播区域,其中K=2π/Λ是空间基波数,Λ是空间基波波长,b(b>1)是尺度参数,h是均方根高度,D(1<D<2)是分形维数,φn是每个谐波的随机初相位,C=是幅值控制参数.图1中展示了锥形波入射无界粗糙表面散射问题的示意图.其中Dr表示由(?)Dr=Hr∪Sr∪r所围区域,Hr是以O为圆心,以r为半径的上半圆弧构成的,Sr={(x,z)|z=f(x),-r<x<r}表示粗糙表面在区间(-r,r)上的限制,Tr由连接Hr与Sr末端的线段组成.另外定义区域Dr,其边界为(?)Dr=Hr∪Sr∪Tr,Hr是以O为圆心,以r为半径的下半圆弧,Tr由连接Hr与Sr末端的线段组成.考虑Thorsos锥形波入射,入射场uinc满足如下非齐次Helmholtz方程其中附加相位项w(r)为且满足其中kg cosθinc》1,θinc(0<|θinc|<π/2)是入射角(波传播方向与z轴负向夹角),g是波束宽度因子,k=2π/λ是入射波波数,λ是入射波波长,i=(?)是虚单位.波束宽度因子g→+∞时,锥形入射波的极限就是平面波.散射场us满足如下Helmholtz方程△us+k2us=0,在Ω内,(3)则全场u满足方程△u+k2u=k2R,在Ω内,(4)及Dirichlet边界条件u|Γ=0.(5)定义0.1设函数us(x,z)∈C2(Ω)满足在Ωl0内,(6)则称函数us满足角谱展开辐射条件,其中uDir(x)∈L1(R)∩L2(R),这里(Fus)(ζ,l0)表示uDir(x)=us(x,l0)(x∈R)的Fourier变换.文中考虑如下锥形波入射无界声软或完全导体表面时散射问题数学模型：问题0.1对于给定的函数f∈C2(R)及由(2)式给定的入射uinc,求满足Helmholtz方程(3)和角谱展开辐射条件(6)的散射场us∈C2(Ω)n C(Ω),使全场u：=us+uinc满足方程(4)和边界条件(5).ⅰ.边界积分方程建立设二维Green函数满足方程其中H0(1)(.)是第一类零阶汉克尔(Hankel)函数,r=(x,z)和r’=(x’,z’)表示R2中的点,δ(·)表示狄拉克(Dirac-δ)函数.引入记号,对于散射问题0.1,为得到与其相应的边界积分方程,首先引入了如下引理：引理0.1若v(r,θ;α)=exp(ikr cos(θ-α)),其中|α|≤π/2,|θ|<π/2,则有I(v；Hr)→0,当r→∞引理0.2若us满足角谱展开辐射条件(6),则有I(us；Hr)→0,当r→∞.引理0.3设平面波up(r,θ;θine):=exp(-ikr cos(θ+θinc)),其中|θinc|<π/2,|θ|<π/2,则有I(up；Hr)→0,I(up；Tr)→0,当r→∞.定理0.1如果全场u满足Helmholtz方程(4)而且散射场us满足角谱展开辐射条件(6),则u满足下面积分表示ⅱ.积分算子的性质引入空间对应范数为|其中定义空间X到空间Y：=PX的算子p：则积分方程能被写成如下算子形式其中定理0.2若f∈C2(R),则算子P：X→Y是有界线性算子对于充分大的L,定义如下形式的算子PT：下面研究算子P和PT的一些性质.定理0.3积分算子PT：X→Y是有界线性算子,并且有Ⅱ.基于快速多极加速的正则化共轭梯度法i.边界积分方程的离散设XT是空间x对应x∈[-L,L]上的限制,同时YT：=PTXT.对应范数为令S：XT→YT为定义0.2定义在{(x,y)|x,y∈(?)D,x≠y}上的K(x,y)称为是弱奇性的,如果K(x,y)在其上连续,且存在M>0和α∈(0,m-1]使得成立,其中(?)D是C1类的有界开区域D(?)Rm的边界.引理0.4C((?)D)是带有模||Φ||∞的Banach空间.定义算子A：C((?)D)→C((?)D),其中核函数K(x,y)是连续的或是弱奇性的,则A在C((?)D)上是紧算子.定理0.4积分算子S：XT→YT是紧线性算子.下面数值求解将积分方程有限截断后所得到的如下形式的积分方程记其算子形式为考虑到Hankel函数的奇性,可得相应积分方程的离散形式：利用Hankel函数的渐进性质,则方程组(12)变为由于算子方程(11)的离散形式中的ρm,在实际数值求解时常常被忽略不计,即求解下述方程组为了得到右端扰动ρm的估计,我们先证明如下的定理.定理0.5方程(13)的右端扰动项有如下估计式成立：其中C是依赖于k,L,θinc,M1,M的常数.ⅱ.快速多极算法的实现在求解线性代数方程组(13)时,配置法的主要限制是计算复杂度因散射面的增长而急剧增加.所以接下来考虑使用快速多极方法,可得相应积分方程的离散形式：从而线性方程组(15)可写成如下矩阵形式其中ⅲ.正则化共轭梯度法尽管研究锥形波入射无界粗糙表面散射问题的相应积分方程时,二重积分是小扰动,但这个小扰动可能对病态问题带来一些困难.更重要的是,它对RCGM的正则化参数选取有重要影响,所以将引进噪声水平δ.采用m=m(δ,bδ)条件终止RCGM,即其中AUmδ-bδ是第m次迭代的残量,bδ是线性方程组(16)的右端扰动项,τ(τ>1)是先验选择常数.对于锥形波入射声软粗糙表面或者完全导体粗糙表面,入射场的能量应该等于粗糙表面的散射能量.所以我们也使用能量关系作为后验迭代停止准则来确保计算的有效性,也就是要求能量指数E=f-π/2π/2dθsσ(θs)=1,其中σ(θs)是双基散射系数.下面给出正则化共轭梯度法.算法0.1.正则化共轭梯度法(RCGM)1.给定噪声水平δ以及最大容许的迭代步数mmax,选取先验常数τ(τ>1),并设m=1；2.给定初始猜测值U0δ,计算r0=bσ-AU0δ,p1=s0=A*r0；3.计算粗糙表面接收到的功率4.设定E0δ=0并给出能量容许误差ε5.只要执行6.否则,终止；6.如果(m≥mmax或||sm||∞≤τδ<||sm-1||∞),终止；否则7.计算8.更新m,m=m+1.值得注意的是算法0.1中的3,4,5,7步是为了保证散射场计算时的能量关系.为了获得停止准则,通过估计ρ(x)得到：由此可以认为,带有余项的积分方程是截断后积分方程的一个扰动.iv.数值结果及其分析本文使用RCGM方法数值求解锥形波入射情形下的一维分形无界粗糙表面的双基散射系数.通过几个具体的数值算例,研究了RCGM的精度和有效性,进一步还考虑了分形粗糙表面的不同参数对散射系数的影响.数值结果表明：(1) RCGM的停止准则对于控制迭代次数是有效的.采用RCGM方法能精确有效的解决锥形波入射情形下大尺度粗糙表面散射问题的计算.这解决了使用矩量法以及经典的共轭梯度法存储量大或缺乏适当停止准则的问题.(2)粗糙表面的参数经常对散射场的分布有非线性的影响,但是我们能看到均方根高度h对双基散射系数的分布有重要的影响,对于给定的入射波,其它参数也对双基散射系数的分布有一些影响.这些结果和经典的Rayleigh准则是一致的.Ⅲ.无界分形粗糙表面多参数反演的迭代法本文中根据散射场信息研究了目标粗糙表面的参数重构和形状反演.这通常对应一个通过探地雷达(GPR)实现地表面(或者海面)成像的模型问题.在GPR系统中,地面上空的发射器向所研究粗糙表面发射电磁波,同时接收器接收数据(见图2),再利用接收到的数据实现表面成像或是目标识别.本文主要研究锥形波入射完全导体(声软)无界粗糙表面参数反演和形状重构.未知参数由最小化目标函数求得,其中目标函数是观测散射场信息和计算散射场差的范数.迭代法中使用RCGM作为正散射问题的快速求解器并结合BFGS优化技术最小化目标函数.ⅰ.反散射问题的快速迭代算法首先考虑已知多重入射波以及相应的散射场信息情况下的粗糙表面参数反演问题的快速迭代算法.第一步考虑M个预先给定的锥形入射波uincm(m=1,2,…,M),照射所要研究的粗糙表面,测量在z=H位置处相应的散射场信息Usm(x,H)(x∈R).当在z=H位置水平安置N个接收器时,在第n个接收器(观测点)处测得的散射场信息为Usm(xn,H)(n=1,2,…,N),简记为Usm,n.第二步v表示粗糙表面中需要重构的未知参数向量.根据粗糙表面未知参数向量的猜测值v0确定相应的粗糙表面.利用本文中的正则化共轭梯度法求解此粗糙表面相应的正散射问题,计算出第n个接收器位置处的散射场信息usm(xn,H)(n=1,2,…,N),简记为usm,n.目标函数表示猜测值为v0时,在所有观测点处散射场的数值解和测量值之差模的平方和的平均值.这里希望找到合适的v使目标函数满足某种误差水平,此时v即为所求.第三步目标函数不为零时,利用迭代算法更新猜测向量v0,最小化目标函数,最终确定向量v0得到待定粗糙表面参数的近似值.这里使用BFGS算法更新猜测向量v0.ⅱ.数值实现和例子目前,由于严格确定粗糙表面参数对散射场的非线性影响仍是有待解决的问题,这导致目标函数关于粗糙表面参数变化的性质缺少理论分析,这是目前迭代法的收敛性难以得到证明的重要原因之一.因此,在研究中通过数值算例检验了逆算法收敛性等性质.同时还进一步讨论了测量误差,初始值,观测点位置以及不同粗糙表面的参数对重构迭代算法的影响.数值实验显示多频多角度入射策略和单一入射策略都能在合理的时间内实现粗糙表面重构.我们可以看到多频多角度入射策略要比单一入射策略时粗糙表面重构的精度高,而且速度相对快.更多还原

【Abstract】 Scattering and inverse scattering of unbounded rough surface is an importantproblem in scientific research and engineering practice, especially in geoscience remotesensing, target recognition, optical difraction and other areas. The theory of scatter-ing and inverse scattering of acoustic and electromagnetic waves plays an importantrole in modern scientific fields. The direct scattering problem, given the informationof the incident wave and the nature of the scatterer, is to find the scattered waveand in particular its behavior at large distances from the scatterer. The inverse scat-tering problem is to find the boundary or physical parameters of the scatterer withgiven incident field and scattered field (or far-field). The scattering problems can besolved numerically by integral equation method, finite element method and infiniteelement method, etc. Roughly speaking, numerical methods for solving inverse scat-tering problems can be classified into three groups: iterative methods, optimization methods and sampling methods. In this dissertation, the study focused on problem of the unbounded rough surface scattering and inverse scattering with tapered wave incident. First, based on the model problem of scattering by unbounded rough surface with tapered wave incident, we strictly deduced the corresponding integral equation and studied the properties of the integral operators. Then, according to the specific characteristics of the corresponding integral equation, we used FMM and RCGM for the numerical solution of the corresponding integral equation. A strategy for the se-lection of the regularization parameter is obtained by estimating the double integral as a perturbed right-hand side of the integral equation. Finally, we discussed the shape reconstruction of the sound soft (perfectly conducting) unbounded rough surface with tapered wave incident. The concrete work is as follows:I. Boundary integral equation methods for scattering by an un-bounded rough surface with tapered wave incidenceLet F:={(x, f(x)) E R2|x E R} denotes the unbounded rough surface, Ω:={r=(x,z)∈R2|z> f(x)} denotes the propagation domain above F, where K=2π/A is the fundamental surface wave number, b(b>1) is the frequency scaling parameter, h is the root mean square height, D(1<D<2) is the fractal dimension, φn is the specified phases, is a normalization factor. The schematic of the unbounded rough surface scattering with tapered wave incidence is shown in Fig.1.The closed region Dr with boundary (?)Dr=Hr U Sr U Tr, where Hr is a large semicircle of radius r and center O in z> f(x), Sr={(x, z)|z=f(x),-r<x<r} is a truncated rough surface, and Tr consists of line segments joining the ends of Hr and Sr. We also define the closed region Dr with boundary (?)Dr=Hr U Sr U Tr, where Hr is a large semicircle of radius r and center O in z<f(x), and Tr consists of the line segments joining the ends of Hr and Sr. Fig.1A schematic of the unbounded rough surface scattering problem with the tapered wave incidenceConsidering the Thorsos tapered wave incident, and the incident field Uinc satisfies nonhomogeneous Helmholtz equation where and additional phase term w(r) is and where kgcosθinc>>1, is the angle of incidence (it is the angle between the direction of propagation and the negative z axis), and g is the parameter that controls the tapering, k=2π/λ is the wave number, λ is the wavelength of incident wave, The tapered incident wave limit is the plane wave as g→+∞. The scattered field us satisfies the Helmholtz equation thus full field u satisfies the equation and Dirichlet boundary conditionsDefinition0.1Let function us∈C2(Ω) is said to satisfy the angular spectrum representation radiation condition (ASRC), if where uDir(x)∈L1(R)∩L2(R), and (Fus)(ζ,l0) is the Fourier transform of uDir(x)=us(x,l0)(x∈R).In this dissertation, we consider the mathematical model of the scattering prob-lem with tapered wave incidence by an unbounded sound soft surface or perfectly conducting surfaces:Problem0.1Give f∈C (R), and an incident field uinc is defined by (2), find us∈C2(Ω)∩C(Ω), such that u:=us+uinc satisfies the equation (4) and boundary condition (5), where us satisfies the Helmholtz equation (3) and the ASRC (6). i. The derivation of boundary integral equationLet the two-dimensional Green’s function satisfies where H0(1)(·) is the first class of zero-order Hankel function, r=(x,z) and r’=(x’,z’) are points in R2,δ(·) is Dirac-δ function. We introduceFor scattering problem0.1, in order to obtain its corresponding boundary integral equation, we firstly prove the following lemmas:Lemma0.1If v(r,θ;α)=exp(ikr cos(θ-α)), where|α|≤π/2, thenLemma0.2If us satisfies the angular spectrum representation radiation con-dition(6), thenLemma0.3Let plane wave up(r,θ;θinc):=exp(-ikr cos(θ+θinc)), where|θinc|<π/2,|θ|<π/2, thenAccording to lemma0.1, lemma0.2and lemma0.3, the corresponding boundary integral equation is derived.Theorem0.1If the total field u satisfies the Helmholtz equation (4), and the scattered field us satisfies the angular spectrum representation radiation condition (6), then u satisfies the integral expression as follows ii. The properties of integral operators We will introduce a space equipped with the corresponding norm where We define the operator P form X to Y:=PX by thus, the integral equation can be written in operator form as whereIn the following, we proved the boundedness of operator P:Theorem0.2For all f∈C2(R), the operator P:X→Y is a bounded linear operator.For sufficiently large L, we define the operator PT byIn the following, we studied the properties of operator P and PT:Theorem0.3The integral operators PT:X→Y is a bounded linear operator, and we haveⅡ. Regularized conjugate gradient method based on the fast multi-pole acceleratedi. Numerical discretization of the boundary integral equationSpace XT is the restriction of X for x on [-L, L], and YT:=PTXT. The relevant norm areLet S:XT→YT as Definition0.2We recall that a kernel K(x,y) is called weakly singular on dD x dD if K(x,y) is denned and continuous for {(x,y)\x,y E dD,x≠y}, and there exist constants M>0and α∈(0,m-1] such that where dD is of class C1.Lemma0.4C((?)D) is a Banach space equipped with the corresponding norm‖·‖∞. Let the mapping A: C(?)D)→C((?)D) be defined by where K(x,y) is continuous on dD x(?)D or weakly singular on dD x dD, Then A defined by (0.4) is compact as an operator from C((?)D) into C((?)D).Theorem0.4Integral operator S: XT→YT is a compact linear operator.Note that the truncated integral equation is that and its operator formApplying the point matching method, and taking into account the singularity of Hankel function and theorem0.1, we get the discretization form of the corresponding integral equation: Applying the small argument approximation of the Hankel function, the equations (12) becomes In the actual numerical solution, for ρm is often ignored in discrete form of operator equation (11), that is solving the following equations: For obtain the estimates of the right-hand member ρm, we proved the following theo-rem.Theorem0.5For the right side of (13), the following estimates hold: where C is a constant which is dependent on k, L,θinc, H1, M1, M.ii. Numerical discretization based on the fast multi-pole acceleratedWhen solving equations (13), the main problem of the above method is that the computational complexity grows dramatically as the size of the scattering surface increases. So we consider to use the FMM, we get the discretization form of the corresponding integral equation: The linear systems (15) can be expressed as whereiii. Regularization conjugate gradient methodCurrently, most of the researches on rough surface scattering problem with tapered wave incidence had omitted the double integral∫Ωk2R(r’)G(r, r’)dσ in the correspond-ing integral equation. Although it is a small perturbation, it is known that a small perturbation for ill-posed problem may bring some troubles. More importantly, it has an impact on the selecting of regularization parameter for RCGM. So we will introduce the noise level S. The RCGM is terminated with m=m(δ, bδ) when where AUmδ—bδ is the residual of the mth iteration, bδ is the perturbed right-hand side of linear equations (16), τ(τ>1) is a priori chosen constant. For scattering of acoustic waves from sound soft surfaces or scattering of horizontally polarized electromagnetic waves from perfectly conducting surfaces, the received power by the rough surface equals to the scattered power. So we also use the energy relation as a posteriori iteration stopping criterion to ensure the effectiveness of computing, i.e. energy index where σ(θs) is the BSC. RCGM is given below.Algorithm0.1.(RCGM)1. Input the noise level S and the largest admissible number of iteration steps mmax, chose constant τ (τ>1), set m=1;2.Input the initial guess U0δ, compute r0=bδ-AU0,δ,p1=s0=A*r0;3.Compute the received power by the rough surface4.Input and the energy stopping tolerance e;5.do while6.if (m≥mmax or‖sm‖∞≤τδ<‖|sm-1‖∞), break; else compute 7. compute8. Set m=m+1;It should be noted that the steps3,4,5,7in algorithm0.1can ensure the energy relation for scattering calculation. In order to obtain the stopping rule, we given the estimate of p(x): Hence, we obtained the noise level.iv. Numerical results and their analysisIn this dissertation, we numerically calculated the BSC for one-dimensional fractal rough surface with tapered wave incidence by using the RCGM for several numerical examples. The influence of different parameters of fractal rough surface on scattering coefficient is considered, the accuracy and efficiency of the RCGM is also studied. Numerical results show that:(1) The stopping criteria of RCGM is effective for controlling the number of it-erations. The RCGM can be performed accurately and efficiently to deal with the scattering problem of large size rough surface with tapered wave incident, which can not be easily accomplished by MOM or classical CGM for their disadvantages of large memory dependence or lack of appropriate stopping criteria.(2) The parameters of rough surface often have nonlinear effects on the distribution of scattering field, we can see that h has significant effect on the distribution of BSC and the other parameters also have some effects on the distribution of BSC for a certain incident wave. In fact, some information of the rough surface itself, such as parameters h, D or K ect., can be retrieved from the data of BSC(or the scattering field). Ⅲ. The iterative method for multi-parameter of unbounded fractal rough surfacesIn this section, we study the parameters reconstruction and shape inversion of the rough surface based on the scattered field. It is a simplification of ground surface (or sea level) imaging via ground-penetrating radar(GPR). In GPR systems, arrays of above-ground transmitters and receivers illuminate surfaces of interest and receive scattered data from the rough surface(see Fig.2), and then using the received data for surface imaging or target recognition. Here, we investigate the parameters inversion and shape reconstruction of unbounded rough surface with tapered wave incidence.Fig.2The schematic of the rough surface inverse scattering with the ta-pered wave incidencei. Fast iterative algorithm of the inverse scattering problemWe presented a new inversion algorithm for the reconstruction of unbounded frac-tal rough surfaces from a set of scattered field measurements for an illumination by tapered wave. The unknown parameters of the surfaces are estimated by minimizing the cost function which is defined as some norm between simulation scattered field data and measured scattered field data. The approach used the fast forward solver-RCGM based on tapered wave incidence combined with the BFGS optimal technique. The flowchart of the inverse scattering algorithm is shown in Fig.3.ii. Numerical implementation and examplesIn general, the cost function is a nonlinear function of the parameters. It is difficult to predict the behavior of cost functions with respect to variations in these parameters. Fig.3Flowchart of inverse scattering algorithmIt is one of the important reasons for the convergence of the iterative method which is difficult to prove. Therefore, in the next study, we examined the nature of the inverse algorithm convergence by numerical examples. The influence of measurement error, initial data, survey position and different rough surfaces on the inversion results were studied. The incidence strategy of multi-angle and multi-frequency can achieve the better results of reconstruction than single incidence within much less iteration steps. The numerical results show that the proposed algorithm is accurate, and the inversion results can be obtained with reasonable computing times.更多还原